!     FILE FOR3D.F                                VERSION 1.1 - 09/04/92
!     ******************************************************************
!     PROGRAM FOR3D
PROGRAM main
!     ******************************************************************
!     *  THIS IS AN IMPROVED VERSION OF THE RESEARCH CODE FOR3D.       *
!     ******************************************************************
!     ******************************************************************
!     *  MATHEMATICAL MODEL BY:                                        *
!     *                                                                *
!     *  DR. D. LEE , CODE 3122                                        *
!     *  ENVIRONMENTAL AND TACTICAL SUPPORT SYSTEMS DEPARTMENT         *
!     *  NAVAL UNDERWATER SYSTEMS CENTER                               *
!     *  NEW LONDON, CT 06320, U.S.A.                                  *
!     *                                                                *
!     *  PROF. M. H. SCHULTZ                                           *
!     *  DEPARTMENT OF COMPUTER SCIENCE                                *
!     *  YALE UNIVERSITY                                               *
!     *  NEW HAVEN, CT 06520, U.S.A.                                   *
!     *                                                                *
!     *  PROF. YOUCEF SAAD                                             *
!     *  DEPARTMENT OF COMPUTER SCIENCE                                *
!     *  UNIVERSITY OF MINNESOTA                                       *
!     *  MINNEAPOLIS, MN 55455-0159, U.S.A.                            *
!     *                                                                *
!     *  PROF. W. L. SIEGMANN                                          *
!     *  DEPARTMENT OF MATHEMATICAL SCIENCES                           *
!     *  RENSSELAER POLYTECHNIC INSTITUTE                              *
!     *  TROY, NY 12180-3590, U.S.A.                                   *
!     ******************************************************************
!     *  ORIGINAL VERSION OF MATHEMATICAL MODEL REPORTED IN:           *
!     *  AN EFFICIENT METHOD FOR SOLVING THE THREE-DIMENSIONAL WIDE    *
!     *  ANGLE WAVE EQUATION, RESEARCH REPORT YALEU/DCS/RR-643,        *
!     *  YALE UNIVERSITY, NEW HAVEN, CT., OCT 1986.                    *
!     ******************************************************************
!     ******************************************************************
!     *  COMPUTER MODEL BY:                                            *
!     *  G. BOTSEAS , CODE 3122                                        *
!     *  ENVIRONMENTAL AND TACTICAL SUPPORT SYSTEMS DEPARTMENT         *
!     *  NAVAL UNDERWATER SYSTEMS CENTER                               *
!     *  NEW LONDON, CONNECTICUT 06320, U.S.A.                         *
!     ******************************************************************
!     *  ORIGINAL VERSION OF COMPUTER MODEL REPORTED IN:               *
!     *  FOR3D: A COMPUTER MODEL FOR SOLVING THE LSS THREE-DIMENSIONAL *
!     *         WIDE ANGLE WAVE EQUATION, NUSC TR 7943, 14AUG87,       *
!     *         G. BOTSEAS, D. LEE, AND D. KING                        *
!     ******************************************************************
!     ******************************************************************
!     *  PROJECT SPONSORED BY:                                         *
!     *  DNL   - PROGRAM MANAGER, DR. G.W. MORTON                      *
!     *  NUSC  - IN-HOUSE LABORATORY INDEPENDENT RESEARCH              *
!     *          PROGRAM MANAGERS, DRS. W.A. VONWINKLE AND K. LIMA     *
!     *  ONR   - PROGRAM MANAGERS, DRS. R.L. LAU, R. OBROCHTA, M. ORR, *
!     *          AND R.BAER                                            *
!     *  NORDA - PROGRAM MANAGER, R.W. MCGIRR                          *
!     ******************************************************************
!     ******************************************************************
!     *  CRAY VERSION - FORTRAN                                        *
!     *  CRAY UNICOS JOB SCRIPT IS IN FILE FOR3D.JOB.                  *
!     ******************************************************************
!     *  MAIN PROGRAM IS FOR3D.                                        *
!     ******************************************************************
!     *  COMMON BLOCK IS IN FILE FOR3D.CMN                             *
!     *  DIMENSIONAL PARAMETERS ARE IN FILE FOR3D.PAR                  *
!ljh  *  Now commons.mod and params.mod
!     ******************************************************************
!     *  ALPHABETICAL LIST OF SUBROUTINES FOLLOWS                      *
!     ******************************************************************
!     *  AMIFD3  - Computes diagonals for matrix A. CALLED BY FOR3D.   *
!     *  BCON3D  - SUPPLIES BOTTOM CONDITION. CALLED BY RHS.           *
!                - IF ITYPEB = 0, BOTX AND BOTY ARE SET TO 0.          *
!                - IF ITYPEB = 1, UBCON3D IS CALLED.                   *
!                - IF ITYPEB = 3, AN ABSORBING LAYER IS INTRODUCED.    *
!     *  BMIFD3  - COMPUTES DIAGONALS FOR MATRIX B. CALLED BY FOR3D.   *
!     *  INDX3D  - COMPUTES INDEX OF REFRACTION ARRAYS XN1 AND XN2.    *
!                - CALLED BY FOR3D.                                    *
!     *  PORT2D  - GENERATES 2D SOLUTION AT PORT SIDEWALL IF ITYPPW=2  *
!     *  PORT3D  - SUPPLIES BOUNDARY CONDITIONS AT PORT SIDEWALL.      *
!                - IF ITYPPW = 0, THE FIELD ALONG THE WALL IS SET TO 0.*
!                - IF ITYPPW = 1, UPORT3D IS CALLED.                   *
!                - IF ITYPPW = 2, PORT3D PREPARES FOR A 2D SOLUTION    *
!                  ALONG THE PORT WALL.                                *
!                - CALLED BY RHS IF NDIM = 3.                          *
!     *  PRINTP  - PRINTS SELECTED INPUT PARAMETERS. CALLED BY FOR3D.  *
!     *  RHS     - COMPUTES RIGHT HAND SIDE OF SYSTEM OF EQUATIONS.    *
!                - CALLED BY FOR3D.                                    *
!     *  SCON3D  - SUPPLIES SURFACE CONDITION. CALLED BY RHS.          *
!                - IF ITYPES = 0, SURX AND SURY ARE SET TO 0.          *
!                - IF ITYPEB = 1, USCON3D IS CALLED.                   *
!     *  SFLD3D  - GENERATES GAUSSIAN STARTING FIELD IF ISF = 0.       *
!                - CALLS USFLD3D IF ISF = 1.                           *
!                - GENERATES GREEN'S STARTING FIELD IF ISF = 2.        *
!                - CALLED BY FOR3D.                                    *
!     *  STBD2D  - GENERATES 2D SOLUTION AT STBD SIDEWALL IF ITYPSW=2  *
!                - AND NDIM = 3. CALLED BY TWOSTEP.                    *
!     *  STBD3D  - SUPPLIES BOUNDARY CONDITIONS AT STARBOARD SIDEWALL. *
!                - IF ITYPSW = 0, THE FIELD ALONG THE WALL IS SET TO 0.*
!                - IF ITYPSW = 1, USTBD3D IS CALLED.                   *
!                - IF ITYPSW = 2, STBD3D PREPARES FOR A 2D SOLUTION.   *
!                  ALONG THE STARBOARD WALL.                           *
!                - CALLED BY RHS IF NDIM = 3.                          *
!     *  SVP3D   - READS SOUND SPEED PROFILES IN INPUT RUNSTREAM.      *
!                - CALLED BY FOR3D IF KSVP = 0.                        *
!     *  TRID3D  - SOLVES A TRIDIAGONAL SYSTEM OF EQUATIONS.           *
!                - CALLED BY PORT2D, STBD2D, AND TWOSTEP.              *
!     *  TWOSTEP - MARCHES THE ACOUSTIC FIELD FORWARD.                 *
!                - FIRST STEP SOLVES IN RANGE.                         *
!                - SECOND STEP SOLVES IN AZIMUTH.                      *
!                - CALLED BY FOR3D.                                    *
!                - IF NDIM = 3 AND ITYPPW = 2, TWOSTEP CALLS PORT2D.   *
!                - IF NDIM = 3 AND ITYPSW = 2, TWOSTEP CALLS STBD2D.   *
!                - IF NDIM = 1 OR 2, TRID3D IS CALLED ONCE.            *
!                - IF NDIM = 3, TRID3D IS CALLED TWICE.                *
!     ******************************************************************
!     *  USER SUBROUTINES                                              *
!     ******************************************************************
!     *  UBCON3D - USER SUPPLIED BOTTOM CONDITIONS. CALLED BY BCON3D.  *
!     *  UBOTTOM - USER SUPPLIED BOTTOM PROFILE INPUT ROUTINE.         *
!                - CALLED BY FOR3D IF KBOT NE 0.                       *
!     *  UEXACT  - EXACT SOLUTION. USER WRITTEN. FOR TEST PURPOSES.    *
!                - CALLED BY FOR3D.                                    *
!     *  UPORT3D - USER SUPPLIED PORT WALL CONDITIONS. CALLED BY PORT3D*
!     *  USCON3D - USER SUPPLIED SURFACE CONDITIONS. CALLED BY SCON3D. *
!     *  USFLD3D - USER SUPPLIED STARTING FIELD.                       *
!                - CALLED BY SFLD3D IF ISF = 1.                        *
!     *  USTBD3D - USER SUPPLIED STBD WALL CONDITIONS. CALLED BY STBD3D*
!     *  USVP3D  - USER SOUND SPEED PROFILE INPUT ROUTINE.             *
!                - CALLED BY FOR3D IF KSVP NOT 0.                      *
!     ******************************************************************
!     *  ALPHABETICAL LIST OF SOME PROGRAM VARIABLES FOLLOWS:          *
!     ******************************************************************
!     *  AL    - ARRAY - Lower diagonal of matrix A - COMPLEX
!     *  ALFA  - SET TO 1/12 IF DOUGLAS METHOD REQUESTED
!     *  ALPHA - VOLUME ATTENUATION - DB/METER
!     *  AM    - ARRAY - Main  diagonal of matrix A - COMPLEX
!     *  ANG   - RELATIVE BEARING - DEGREES
!     *  AU    - ARRAY - Upper diagonal of matrix A - COMPLEX
!     *  ATT   - Attenuation coefficient for artificial absorbing layer
!     *  BETA  - ARRAY - ATTENUATION IN LAYERS - DB/WAVELENGTH
!     *  BL    - ARRAY - LOWER DIAGONAL OF MATRIX B - COMPLEX
!     *  BM    - ARRAY - MAIN DIAGONAL OF MATRIX B - COMPLEX
!     *  BOTX  - ARRAY - COMPLEX PRESSURE AT BOTTOM AT ADVANCED RANGE RA+DR
!     *  BOTY  - ARRAY - COMPLEX PRESSURE AT BOTTOM AT PRESENT RANGE RA
!     *  BTA   - ARRAY - PARTIAL SOLUTION OF SYSTEM OF EQUATIONS
!     *  BU    - ARRAY - UPPER DIAGONAL OF MATRIX B - COMPLEX
!     *  C0    - REFERENCE SPEED OF SOUND - METERS/SEC
!     *  CSVP  - ARRAY - SOUND VELOCITY - METERS/SEC
!     *  D     - ARRAY - RIGHT HAND SIDE OF SYSTEM OF EQUATIONS - COMPLEX
!     *  DEG   - CONVERSION FACTOR - DEGREES/RADIAN
!     *  DOUGRA - 0     - DON'T USE DOUGLAS METHOD
!                 NOT 0 - RANGE AT WHICH TO BEGIN DOUGLAS METHOD - METERS
!     *  DR    - RANGE STEP - METERS
!     *  DTH   - WIDTH OF SECTOR - DEGREES
!     *  DZ    - DEPTH INCREMENT OF SOLUTION - METERS
!     *  DZZ   - DEPTH INCREMENT FOR ADJUSTING LAYER DEPTHS IN SLOPING BOTTOM
!     *  FLDW  - WIDTH OF FIELD - DEGREES
!     *  FRQ   - FREQUENCY - HZ
!     *  HNK   - HANKEL FUNCTION H0(1)
!     *  HNKL  - EXTERNAL FUNCTION - COMPUTES HANKEL FUNCTION H0(1)
!     *  IBOT  - BOTTOM DEPTH PRINT FLAG
!                0 - DO NOT PRINT BOTTOM DEPTHS
!                1 - PRINT BOTTOM DEPTHS
!     *  IHNK  - HANKEL FUNCTION FLAG
!                0 - HANKEL FUNCTION NOT USED. 10*LOG(R) ADDED TO SOLUTION.
!                1 - STARTING FIELD DIVIDED BY HANKEL FUNCTION.
!                    SOLUTION MULTIPLIED BY HANKEL FUNCTION BEFORE
!                    COMPUTING PROPAGATION LOSS.
!     *  ILYR  - INDEX FOR ARRAYS BETA, ZLYR, AND RHO
!     *  IPZ   - EVERY IPZ'TH VALUE IN DEPTH IS PRINTED
!     *  ISF   - STARTING FIELD FLAG
!                0 - PROGRAM GENERATES GAUSSIAN STARTING FIELD
!                    AT RANGE = 0.0. SEE SUBROUTINE SFLD3D.
!                1 - USER SUPPLIES STARTING FIELD. SEE SUBROUTINE USFLD3D.
!                2 - GREENS WIDE ANGLE STARTER
!                3 - AVAILABLE
!     *  ISFLD - STARTING FIELD PRINT FLAG
!                0 - DO NOT PRINT STARTING FIELD
!                1 - PRINT STARTING FIELD
!     *  ISVP  - SVP PRINT FLAG
!                0 - DO NOT PRINT SOUND VELOCITY PROFILE
!                1 - PRINT SOUND VELOCITY PROFILE
!     *  ITEMP - TEMPORARY VARIABLE
!     *  ITRK  - ARRAY - INDICES FOR ARRAY TRACK
!     *  ITYPEB - TYPE OF BOTTOM
!                 0 - BCON3D SETS BOTY AND BOTX = 0.0
!                 1 - USER SUPPLIES BOTTOM CONDITION. SEE SUBROUTINE UBCON3D.
!                 2 - SPARE.
!                 3 - ABSORBING LAYER INTRODUCED - BOTTOM OF LAYER IS FLAT
!                 4 - SPARE.
!     *  ITYPES - TYPE OF SURFACE
!                 0 - PRESSURE RELEASE. SCON3D SETS SURY AND SURX = 0.0
!                 1 - USER SUPPLIES SURFACE CONDITION. SEE SUBROUTINE USCON3D.
!                 2 - SPARE
!     *  ITYPPW - TYPE OF PORT SIDEWALL BOUNDARY CONDITION.
!                 0 - FIELD ALONG PORT SIDEWALL IS SET TO 0.0.
!                 1 - USER SUPPLIED. SEE SUBROUTINE UPORT3D.
!                 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3.
!     *  ITYPSW - TYPE OF STBD SIDEWALL BOUNDARY CONDITION.
!                 0 - FIELD ALONG STBD SIDEWALL IS SET TO 0.0.
!                 1 - USER SUPPLIED. SEE SUBROUTINE USTBD3D.
!                 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3.
!     *  IWZ    - INDEX INCREMENT OF RECEIVER SOLUTIONS TO BE WRITTEN ON DISK.
!     *  IXSVP  - ARRAY OF POINTERS WHICH POINT TO ENTRIES IN CSVP AND ZSVP.
!               - IXSVP(1,L) POINTS TO BOTTOM DEPTH AND SPEED IN LAYER 1
!                 ON SECTOR BOUNDARY L.
!               - IXSVP(2,L) POINTS TO BOTTOM DEPTH AND SPEED IN LAYER 2
!                 ON SECTOR BOUNDARY L.
!                 ETC.
!     *  KSVP   - SVP PROFILE FLAG.
!                 0      - PROFILE IS IN RUNSTREAM - SEE SUBROUTINE SVP3D.
!                 NOT 0  - USER ROUTINE GENERATES PROFILES - SEE SUBROUTINE
!                          USVP3D.
!               - KSVP MAY BE USED AS AN INDEX IN A GOTO STATEMENT IN USVP3D.
!     *  KBOT   - BOTTOM PROFILE FLAG.
!                 0      - PROFILE IS IN RUNSTREAM .
!                 NOT 0  - USER ROUTINE GENERATES PROFILES - SEE SUBROUTINE
!                          UBOTTOM.
!     *  MLYR  - TEMPORARY - NUMBER OF LAYERS INVOLVED IN SPECIFIC CALCULATION.
!     *  MM    - INDEX - MM+1 POINTS TO FIRST VALUE OF ARTIFICIAL ABSORBING
!                        LAYER IN ARRAY U.
!     *  MXLYR - PARAMETER - MAXIMUM NUMBER OF LAYERS.
!                           - MAX DIMENSION OF ARRAYS BETA,RHO,ZLYR,IXSVP.
!     *  MXN   - PARAMETER - MAXIMUM DIMENSION OF C, X, Y, AND U ARRAYS.
!     *  MXSVP - PARAMETER - MAXIMUM DIMENSION OF ARRAYS CSVP AND ZSVP.
!     *  MXTRK - PARAMETER - MAXIMUM DIMENSION OF ARRAY TRACK.
!     *  N     - NUMBER OF EQUI-SPACED POINTS IN U.
!              - INCLUDES BOTTOM POINT - DOES NOT INCLUDE SURFACE POINT.
!     *  N1    - Diagonal elements N1 thru N will be computed - not used.
!     *  NA    - NUMBER OF POINTS IN ABSORBING LAYER.
!     *  NDIV  - IF DOUGLAS METHOD IS REQUESTED, DIVIDE N BY NDIV.
!     *  NIU   - PARAMETER - INPUT UNIT NUMBER.
!     *  NLYR  - NUMBER OF LAYERS.
!     *  NOLD  - NUMBER OF RECEIVER DEPTHS ON ENTRY TO ROUTINE CRNK.
!     *  NOU   - PARAMETER - NPU_ENABLE UNIT NUMBER.
!     *  NPU   - PARAMETER - PRINTER UNIT NUMBER.
!     *  NSEC  - NUMBER OF SECTORS.
!     *  NSOL  - NUMBER OF SOLUTIONS (NSEC+1).
!     *  NSVP  - NUMBER OF POINTS IN CSVP AND ZSVP.
!     *  NTOT  - TOTAL POINTS IN SOLUTION (N*NSOL).
!     *  NWSVP - NUMBER OF POINTS IN LAYER 1 SVP.
!     *  NZ    - NUMBER OF SOLUTION DEPTHS TO BE WRITTEN ON DISK.
!     *  OLDR  - Range increment at start of problem.
!     *  PDR   - RANGE INCREMENT AT WHICH SOLUTION IS PRINTED - METERS.
!     *  PDTH  - AZIMUTHAL INCREMENT AT WHICH SOLUTION IS PRINTED - DEGREES.
!     *  PDZ   - DEPTH INCREMENT AT WHICH SOLUTION IS PRINTED - METERS.
!     *  PHI   - Width of sector in radians.
!     *  PI    - 3.141592654
!     *  PL    - PROPAGATION LOSS - DB.
!     *  PORTX - ARRAY - COMPLEX FIELD AT PORT SIDEWALL - ADVANCED RANGE.
!     *  PORTY - ARRAY - COMPLEX FIELD AT PORT SIDEWALL - PRESENT RANGE.
!     *  R0    - INITIAL RANGE - METERS.
!     *  R1    -      RANGE AT WHICH BOTTOM DEPTH IS AVAILABLE - METERS.
!     *  R2    - NEXT RANGE AT WHICH BOTTOM DEPTH IS AVAILABLE - METERS.
!     *  RA    - HORIZONTAL RANGE OF STARTING FIELD FROM SOURCE - METERS.
!              - RA IS SET TO 0.0 IF STARTING FIELD IS GAUSSIAN. RA IS
!              - INCREMENTED BY DR AS SOLUTION IS MARCHED OUT IN RANGE.
!     *  RA+DR - RANGE TO WHICH SOLUTION IS TO BE ADVANCED - METERS.
!     *  RHO   - ARRAY - Density in layer.
!     *  RHOG  - ARRAY - Density gradient in layer.
!     *  RMAX  - MAXIMUM RANGE OF SOLUTION - METERS.
!     *  RSVP  - NEXT RANGE AT WHICH SVP IS AVAILABLE - METERS.
!     *  STBDX - ARRAY - COMPLEX FIELD AT STBD SIDEWALL - ADVANCED RANGE.
!     *  STBDY - ARRAY - COMPLEX FIELD AT STBD SIDEWALL - PRESENT RANGE.
!     *  SURX  - ARRAY - COMPLEX PRESSURE AT SURFACE AT ADVANCED RANGE RA+DR.
!     *  SURY  - ARRAY - COMPLEX PRESSURE AT SURFACE AT PRESENT RANGE RA.
!     *  TEMP  - TEMPORARY VARIABLE - COMPLEX.
!     *  THETA - ARRAY - SLOPE OF BOTTOM ALONG SECTOR BOUNDARIES - RADIANS.
!                .EQ.0 -- FLAT BOTTOM .
!                .GT.0 -- SHALLOW TO DEEP.
!                .LT.0 -- DEEP TO SHALLOW.
!     *  TITLE - 80 CHARACTER TITLE.
!     *  TM    - ARRAY - TIME OF DAY.
!     *  TRACK - 3 DIM. ARRAY - RANGE AND DEPTH OF WATER - METERS.
!              - ASSUMES FLAT BOTTOM AT PRESENT TIME.
!     *  U     - ARRAY - COMPLEX ACOUSTIC PRESSURE FIELD.
!     *  U1-U12 - USER MAY USE THESE VARIABLES TO INPUT DATA REQUIRED BY USER
!               - ROUTINES - REAL - IN COMMON BLOCK.
!     *  UA-UE  - USER MAY USE THESE VARIABLES - COMPLEX - IN COMMON BLOCK
!     *  V      - ARRAY - PARTIAL SOLUTION OF COMPLEX FIELD
!     *  WDR    - RANGE STEP AT WHICH SOLUTION IS WRITTEN ON DISK - METERS
!     *  WDTH   - AZIMUTHAL INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK - DEG
!     *  WDZ    - DEPTH INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK - METERS
!               - WDZ SHOULD BE SELECTED SO THAT PLOT PROGRAM DOES NOT
!               - INTERPOLATE BETWEEN WIDELY SPACED RECEIVERS.
!     *  WZ1    - FIRST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK
!     *  WZ2    - LAST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK
!               - IN OTHER WORDS, WRITE WZ1 TO WZ2 BY WDZ - METERS
!     *  XK0    - Reference Wave Number
!     *  XN2    - ARRAY - INDEX OF REFRACTION SQUARED - COMPLEX
!     *  XPR    - RANGE AT WHICH SOLUTION IS PRINTED - METERS
!     *  XWR    - RANGE AT WHICH SOLUTION IS WRITTEN ON DISK - METERS
!     *  Z1     - DEPTH OF WATER AT RANGE R1 - METERS
!     *  Z2     - DEPTH OF WATER AT RANGE R2 - METERS
!     *  ZA     - DEPTH OF FIELD AT RANGE RA - METERS
!               - INITIAL DEPTH OF STARTING FIELD AT RANGE RA IS
!               - AS FOLLOWS:
!               - IF ITYPEB = 0 OR 1, ZA IS MAXIMUM DEPTH OF
!               - BOTTOM-MOST SEDIMENT LAYER AT INITIAL RANGE OF
!               - STARTING FIELD. IF ITYPEB = 3, ZA IS MAXIMUM
!               - DEPTH OF ARTIFICIAL ABSORBING LAYER AT INITIAL
!               - RANGE OF STARTING FIELD. PROGRAM INSERTS LAYER.
!               - RHO AND BETA ARE OBTAINED FROM LAYER ABOVE.
!               - SPEED IS BOTTOM-MOST SPEED FROM LAYER ABOVE.
!               - BOTTOM OF ABSORBING LAYER REMAINS FLAT.
!     *  ZI     - DEPTH OF RECEIVER 'I' - METERS.
!     *  ZLYR   - ARRAY - DEPTH TO BOTTOM OF LAYER - METERS.
!     *  ZS     - SOURCE DEPTH - METERS.
!     *  ZSVP   - ARRAY - DEPTH OF SOUND VELOCITY - METERS.
!     ******************************************************************
!     ******************************************************************
!     *  INPUT UNIT NUMBER = NIU
!     *  INPUT FILE NAME   = for3d.in - SEE CRAY UNICOS JOB SCRIPT
!     *  CONTENTS: INPUT RUNSTREAM IN FREE FORMAT
!
!  LINE    INPUT DATA
!
!    1     TITLE
!
!    2     NDIM
!
!    3     FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,ITYPPW,ITYPSW,FLDW,NSEC
!
!    4     RMAX,DR,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH,ISFLD,ISVP,IBOT
!
!    5     DOUGRA,NDIV
!
!    6     U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
!
!    7     KBOT
!
!    8     R1,Z1  **     IF KBOT = 0, BOTTOM PROFILE IS IN RUNSTREAM.
!    9     R2,Z2    *    RANGE, WATER DEPTH (METERS).
!            .      **   -1,-1 MARKS END OF THIS PROFILE.
!            .      ***  ENTER PROFILE FOR EACH SECTOR BOUNDARY.
!            .      **   PORT BOUNDARY FIRST.
!    N       .      *    IF KBOT NE 0, UBOTTOM IS CALLED. OMIT
!    N+1   -1,-1  **     LINES 8 THRU N+1.
!
!             ************************************************
!    N+2  RSVP                                                *
!    N+3  KSVP                                                *
!             **************************************          *
!    N+4  NLYRS(L)                                  *         *
!             ****************************          *         **
!    N+5  ZLYR(I,L),RHO(I,L),RHOG(I,L),   *         **        *** NOTE
!                   BETA(I,L),BETAG(I,L)  **        *** NOTE  ***  1
!    N+6  ZSVP(1,L),CSVP(1,L)             *** NOTE  ***  2    **
!    N+7  ZSVP(2,L),CSVP(2,L)             **   3    **        *
!     .          .       .                *         *         *
!    N+M  ZSVP(J,L),CSVP(J,L) ************          *         *
!                             **********************          *
!                             ********************************
!
!     ******************************************************************
!     *** QUICK REFERENCE AND NOTES FOR INPUT RUNSTREAM
!     *** UNITS: METERS AND METERS/SEC EXCEPT AS NOTED
!     ******************************************************************
!  LINE   INPUT DATA
!
!   1     TITLE  =  USER NOTE. 80 CHARACTERS MAXIMUM.
!
!   2     NDIM   =  NUMBER OF DIMENSIONS.
!                   1 - 2D SOLUTION.
!                   2 - N x 2D SOLUTIONS.
!                   3 - 3D SOLUTIONS.
!
!   3     FRQ    =  FREQUENCY (HZ)
!
!         ZS     =  SOURCE DEPTH
!
!         C0     =  REFERENCE SOUND SPEED. IF C0 = 0.0, C0 IS SET TO AVERAGE
!                   SPEED IN FIRST LAYER.
!
!         ISF    =  STARTING FIELD FLAG.
!                   0 = GAUSSIAN.
!                   1 = USER FIELD.
!                   2 = GREENS STARTER.
!
!         RA     =  HORIZONTAL RANGE FROM SOURCE TO STARTING FIELD.
!                   RA IS SET TO 0.0 IF ISF = 0.
!
!         ZA     =  DEPTH OF STARTING FIELD AT RANGE RA. IF ZA = 0.0, ZA IS SET
!                   TO MAX DEPTH OF BOTTOM LAYER IN FIRST PROFILE. IF ITYPEB =
!                   2 OR 3 AND ZA = 0.0, ZA IS SET TO (4/3)*MAX DEPTH OF BOTTOM
!                   LAYER. IF ITYPEB = 3 AND ZA NOT ZERO, THE ARTIFICIAL
!                   BOTTOM LAYER IS EXTENDED TO ZA METERS PROVIDED THAT ZA
!                   IS GREATER THAN OR EQUAL TO MAX DEPTH OF BOTTOM LAYER
!                   IN FIRST PROFILE.
!
!         N      =  NUMBER OF EQUISPACED RECEIVERS IN STARTING FIELD. IF N = 0,
!                   N IS SET SO THAT THE RECEIVER DEPTH INCREMENT IS
!                   EQUAL TO 1 METER. IF N IS GREATER THAN MXN
!                   N IS SET TO MXN.
!
!         IHNK   =  HANKEL FUNCTION FLAG. IHNK = 0, DON'T USE HANKEL FUNCTION.
!                   IHNK = 1, DIVIDE STARTING FIELD BY HANKEL FUNCTION, THEN
!                   MULTIPLY THE SOLUTION FIELD BY HANKEL FUNCTION BEFORE
!                   COMPUTING PROPAGATION LOSS. IF STARTING FIELD IS GAUSSIAN,
!                   IHNK SHOULD BE SET TO 0. IF STARTING FIELD IS ELLIPTIC,
!                   IHNK SHOULD BE SET TO 1.
!
!         ITYPES =  TYPE OF SURFACE
!                   0 - PRESSURE RELEASE. SCON3D SETS SURY AND SURX = 0.0
!                   1 - USER SUPPLIES SURFACE CONDITION. SEE SUBROUTINE USCON3D.
!                   2 - SPARE.
!
!         ITYPEB =  TYPE OF BOTTOM
!                   0 - PRESSURE RELEASE. BCON3D SETS BOTY AND BOTX = 0.0
!                   1 - USER SUPPLIES BOTTOM CONDITION. SEE SUBROUTINE UBCON3D.
!                   2 - SPARE.
!                   3 - ABSORBING LAYER INTRODUCED - FLAT BOTTOM
!                   4 - SPARE
!
!         ITYPPW =  TYPE OF PORT SIDEWALL BOUNDARY CONDITION
!                   0 - FIELD ALONG PORT SIDEWALL IS SET TO 0.0
!                   1 - USER SUPPLIED. SEE SUBROUTINE UPORT3D.
!                   2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3.
!
!         ITYPSW =  TYPE OF STARBOARD SIDEWALL BOUNDARY CONDITION
!                   0 - FIELD ALONG STARBOARD SIDEWALL IS SET TO 0.0
!                   1 - USER SUPPLIED. SEE SUBROUTINE USTBD3D.
!                   2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3.
!
!         FLDW   =  WIDTH OF FIELD IN DEGREES. IGNORED IF NDIM = 1.
!
!         NSEC   =  NUMBER OF SECTORS IN FIELD. IGNORED IF NDIM = 1.
!                   NUMBER OF SOLUTIONS = NSOL = NSEC+1.
!
!  4      RMAX   =  MAXIMUM RANGE OF SOLUTION
!
!         DR     =  RANGE STEP. IF DR = 0, DR IS SET TO 1 METER.
!                   IF BOTTOM OF PROBLEM IS NOT FLAT, DR IS RECOMPUTED
!                   SO THAT MAX DEPTH IS EITHER INCREMENTED OR DECREMENTED
!                   BY DZ. SOLUTION IS COMPUTED EVERY DR METERS.
!
!         WDR    =  RANGE STEP AT WHICH SOLUTION IS WRITTEN ON DISK.
!                   IF WDR NOT 0, AN NPU_ENABLE DISK FILE IS ASSIGNED.
!                   WDR IS ROUNDED TO NEAREST DR.
!
!         WZ1    =  FIRST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK.
!
!         WZ2    =  LAST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK.
!
!         WDZ    =  DEPTH INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK.
!                   ROUNDED TO NEAREST DZ.
!
!         WDTH   =  AZIMUTHAL INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK.
!                   ROUNDED TO NEAREST DTH.
!
!         PDR    =  RANGE STEP AT WHICH SOLUTION IS PRINTED.
!                   ROUNDED TO NEAREST DR.
!
!         PDZ    =  DEPTH INCREMENT AT WHICH SOLUTION IS PRINTED.
!                   ROUNDED TO NEAREST DZ.
!
!         PDTH   =  AZIMUTHAL INCREMENT AT WHICH SOLUTION IS PRINTED.
!                   ROUNDED TO NEAREST DTH.
!
!         ISFLD  =  0 - DON'T PRINT STARTING FIELD.
!                =  1 - PRINT STARTING FIELD.
!
!         ISVP   =  0 - DON'T PRINT SOUND VELOCITY PROFILE.
!                =  1 - PRINT SOUND VELOCITY PROFILE.
!
!         IBOT   =  0 - DON'T PRINT BOTTOM DEPTHS.
!                =  1 - PRINT BOTTOM DEPTHS.
!
!  5      DOUGRA =  RANGE AT WHICH TO SWITCH FROM CRANK-NICOLSON METHOD TO
!                   DOUGLAS METHOD. IF DOUGRA = 0, USE CRANK-NICOLSON METHOD
!                   ONLY. SUGGESTED VALUE IS 5000 METERS.
!
!         NDIV   =  IF DOUGRA NE 0, N IS DIVIDED BY NDIV RESULTING IN ONLY
!                   N/NDIV SOLUTIONS IN DEPTH. SPEED UP IS NDIV TIMES.
!                   SUGGESTED VALUE IS 5.
!
!  6      U1-U12 =  USER VARIABLES - REAL, SINGLE PRECISION.
!                   SEE HARVARD SUBROUTINE HARVARD.F FOR EXAMPLE.
!
!  7      KBOT   =  BOTTOM PROFILE FLAG.
!                =  0 - BOTTOM PROFILE IN INPUT RUNSTREAM.
!                =  NOT ZERO. PROFILE IS SUPPLIED BY SUBROUTINE UBOTTOM.
!                   USER WRITES UBOTTOM. IF NSOL IS LARGE, UBOTTOM IS
!                   PREFERRED METHOD OF INPUT.
!
!                  BOTTOM PROFILE AT LEFTMOST SECTOR BOUNDARY.
!  8      R1,Z1  = RANGE AND DEPTH OF WATER.
!  9      R2,Z2  = ETC.
!           .
!  N+1    -1,-1  = MARKS THE END OF THIS BOTTOM PROFILE.
!
!                  NEXT PROFILE.
!  8      R1,Z1  = RANGE AND DEPTH OF WATER.
!  9      R2,Z2  = ETC.
!           .
!  N+1    -1,-1  = MARKS THE END OF THIS BOTTOM PROFILE.
!
!                  ETC. REPEAT UNTIL ALL PROFILES ENTERED.
!
!  N+2    RSVP   = RANGE OF SOUND SPEED PROFILES.
!                  SEE NOTE 1.
!
!  N+3    KSVP   = SVP FLAG.
!                = 0 - SOUND SPEED PROFILES IN INPUT RUNSTREAM.
!                = NOT ZERO.  PROFILE (LINES N+4 THRU N+M) IS SUPPLIED BY
!                  USER. USER WRITES SUBROUTINE USVP3D. KSVP MAY BE USED IN
!                  COMPUTED GOTO STATEMENT TO TRANSFER CONTROL IN USVP3D.
!                  IF NSOL IS LARGE, USVP3D IS PREFERRED METHOD OF INPUT.
!
!  N+4    NLYRS(L) = NUMBER OF LAYERS. IF ITYPEB  = 3, PROGRAM INSERTS
!                    AN ARTIFICIAL LAYER AND INCREMENTS NLYRS(L) BY 1.
!                    SEE NOTE 2.
!
!  N+5    ZLYR(I,L) = MAX DEPTH OF LAYER I IN PROFILE.
!
!         RHO(I,L)  = DENSITY IN LAYER I (G/CM**3).
!
!         RHOG(I,L) = DENSITY GRADIENT IN LAYER I (G/CM**3/M).
!
!         BETA(I,L) = ATTENUATION IN LAYER I (DB/WAVELENGTH).
!                     IF BETA(I,L) IS NEGATIVE, ATTENUATION IS COMPUTED.
!                     SEE NOTE 3.
!
!         BETAG(I,L)= ATTENUATION GRADIENT IN LAYER I (DB/WAVELENGTH/M).
!
!  N+6    ZSVP(1,L) = DEPTH TO TOP OF LAYER I
!
!         CSVP(1,L) = SPEED OF SOUND AT TOP OF LAYER I
!            .
!            .
!  N+M    ZSVP(J,L) = DEPTH TO BOTTOM OF LAYER I
!
!         CSVP(J,L) = SPEED OF SOUND AT BOTTOM OF LAYER I
!                     IF ONLY ONE SVP INPUTTED, IT IS USED THRU ENTIRE PROBLEM.
!                     IF MORE THAN ONE SVP INPUTTED, LAST SVP IS USED THRU
!                     REMAINDER OF PROBLEM.
!
!  NOTE 1  REPEAT LINES N+2 THRU N+M FOR EACH SET OF PROFILES TO BE ENTERED.
!          A SET OF PROFILES CONSISTS OF NSOL (NSEC+1) PROFILES LOCATED ON
!          SECTOR BOUNDARIES AND EQUISPACED FROM PORT TO STARBOARD IN AZIMUTH.
!          NSEC IS THE NUMBER OF SECTORS IN THE FIELD OF INTEREST. NSEC
!          ADJACENT SECTORS HAVE NSOL BOUNDARIES.
!
!  NOTE 2  REPEAT LINES N+4 THRU N+M FOR EACH PROFILE IN THE SET. AT RANGE ZERO
!          ALL PROFILES ARE THE SAME. BEYOND RANGE ZERO, L PROFILES ARE
!          ENTERED FROM PORT TO STARBOARD AS VIEWED FROM RANGE R0. L IS
!          VARIED FROM 1 THRU NSOL. ENTERING -1,-1 AFTER PROFILE L
!          WILL CAUSE THAT PROFILE TO BE ENTERED REPEATEDLY UNTIL
!          NSOL PROFILES HAVE BEEN ENTERED.
!
! NOTE 3  REPEAT LINES N+5 THRU N+M FOR EACH LAYER I IN THE PROFILE.
!         J IS THE NUMBER OF POINTS IN LAYER I.
!
!     ******************************************************************
!     *** NPU_ENABLE
!     ******************************************************************
!         NPU_ENABLE UNIT NUMBER = NOU
!         NPU_ENABLE FILE NAME   = for3d.out - SEE CRAY UNICOS JOB SCRIPT
!         CONTENTS:  AS FOLLOWS - UNFORMATTED
!
!         THE FOLLOWING VARIABLES ARE WRITTEN ONCE:
!
!   WRITE(NOU)NDIM,FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,ITYPPW,ITYPSW,FLDW,
!             NSEC,NSOL,RMAX,DR,WDR,WZ1,WZ2,WDZ,DZ,DOUGRA,NDIV,
!             U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
!
!    THE VALUES OF THE VARIABLES AND ARRAYS LISTED IN THE FOLLOWING
!    WRITE STATEMENT ARE WRITTEN AT R0 AND AT EACH WDR (WRITE RANGE INCREMENT)
!    THEREAFTER.
!
!   WRITE(NOU)ANG,NZ,RA,WZ1,WDZ,(U(M+I),I=IWZ1,IWZ2,IWZ)
!
!   where  M   = (J-1)*N
!          J   = SOLUTION INDEX. J VARIES FROM 1 TO NSOL. PORT TO STARBOARD.
!          N   = NUMBER OF RECEIVERS IN VERTICAL COLUMN.
!          ANG = RELATIVE BEARING OF SOLUTION IN DEGREES. PORT IS NEGATIVE.
!          NZ  = NUMBER OF RECEIVERS WRITTEN ON DISK.
!          RA  = RANGE IN METERS.
!          WZ1 = DEPTH OF FIRST RECEIVER.
!          WDZ = DEPTH INCREMENT OF RECEIVERS WRITTEN ON DISK IN METERS.
!          U   = SOLUTION ARRAY. COMPLEX.
!
!     PRINTED NPU_ENABLE IS DIRECTED TO PRINTER UNIT NUMBER NPU.
!     ******************************************************************
!     ******************************************************************
!     ******************************************************************
!
!     *** INCLUDE COMMON BLOCK
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      IMPLICIT NONE
      CHARACTER(LEN=80) :: TITLE               ! cljh
      COMPLEX           :: HNKL            ! cljh
      REAL    (KIND=8)  :: ATT, DRLAST, DZZ, &
                           OLDR, OLDZ, PL, RTEMP, &
                           TH, XPR, XWR, ZBOT, ZG, ZMAX, ZWSVP
      INTEGER (KIND=4)  :: I, ILYR, IPZ, &
                           ITEMP, IWZ, IWZ1, IWZ2, &
                           J, JSOL, K, KBOT, L, M, MLYR, MM, &
                           MOLD, NA, NOLD, NSOLP, &
                           NSOLW, NTOLD, &
                           NTRK, NZ
!
!     *** READ INPUT PARAMETERS.
!     *** FILES ARE LINKED TO UNIT NUMBERS IN CRAY JOB SCRIPT.
      OPEN(UNIT=NIU,FILE='FOR3D.IN')
!
      READ (NIU, 15 ) TITLE
      WRITE (NPU, 15 ) TITLE
   15 FORMAT(A80)
!
!     *** READ SOLUTION DIMENSION
      READ(NIU,*)NDIM
      IF (NDIM < 1 .OR. NDIM > 3) THEN
        WRITE(NPU,*)'INPUT ERROR. CHECK NDIM VALUE.'
        STOP
      ELSE
      ENDIF
!
      READ(NIU,*)FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,    &
      ITYPPW,ITYPSW,FLDW,NSEC
!
      READ(NIU,*) RMAX,DR,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH,    &
      ISFLD,ISVP,IBOT
!
      READ(NIU,*) DOUGRA,NDIV
!
      READ(NIU,*) U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
!
!     SET DOUGLAS COEFFICIENT TO 0
      ALFA=0.0
!
!     *** IF C0 NE 0.0 COMPUTE REFERENCE WAVE NUMBER
      IF (C0 /= 0.0)XK0=2.0*PI*FRQ/C0
!
!     *** COMPUTE NUMBER OF SOLUTIONS. ONE FOR EACH SECTOR BOUNDARY.
!
      IF (NDIM == 1) THEN
!       *** USER REQUESTED A 2D SOLUTION. ONLY ONE SOLUTION IS GENERATED.
        FLDW=0.0
        NSEC=0
        NSOL=1
        NSOLW=1
        NSOLP=1
        DTH=0
        WDTH=0.0
        PDTH=0.0
        ANG=0.0
        PHI=0.0
!
      ELSE
!       *** MORE THAN 1 SOLUTION.
        NSOL=NSEC+1
!       *** COMPUTE SECTOR WIDTH IN DEGREES.
        DTH=FLDW/NSEC
!       *** COMPUTE SECTOR WIDTH IN RADIANS.
        PHI=PI*FLDW/(NSEC*180)
!       *** COMPUTE NUMBER OF SOLUTIONS IN AZIMUTH TO WRITE TO DISK.
!       *** AVOID CRAY ROUND UP PROBLEMS
        NSOLW=INT(WDTH/DTH+.5)
        IF (NSOLW == 0)NSOLW=1
!       *** COMPUTE NUMBER OF SOLUTIONS IN AZIMUTH TO PRINT.
        NSOLP=INT(PDTH/DTH+.5)
        IF (NSOLP == 0)NSOLP=1
      ENDIF
!
      IF (NSOL > MXSOL) THEN
        WRITE(NPU,*)'TOO MANY SOLUTIONS REQUESTED.'
        WRITE(NPU,*)'INCREASE PARAMETER MXSEC AND RECOMPILE.'
        STOP
      ELSE
      ENDIF
!
!     *** IF GAUSSIAN STARTING FIELD, RA MUST BE 0.
      IF (ISF == 0) RA=0.0
      R0=RA
!
!     *** READ BOTTOM PROFILE FLAG
      READ(NIU,*) KBOT
!
!     *** IS PROFILE IN RUNSTREAM?
      IF (KBOT == 0) THEN
!       *** YES
!       *** READ BOTTOM PROFILE - RANGE,DEPTH
        ZBOT=0.0
        DO 24 L=1,NSOL
          DO 22 I=1,MXTRK+1
            READ(NIU,*) TRACK(I,1,L),TRACK(I,2,L)
!           *** END OF PROFILE?
            IF (TRACK(I,1,L) < 0.0) GO TO 23
!           *** NO
            NTRK=I
            IF (ZBOT < TRACK(I,2,L))ZBOT=TRACK(I,2,L)
   22     CONTINUE
   23     CONTINUE
!
!         *** ERROR DETECTED?
          IF (NTRK == 1 .OR. NTRK > MXTRK) THEN
!           *** YES.
            ITEMP=MXTRK
            WRITE(NPU,29) ITEMP
            STOP
          ELSE
!           *** NO ERROR DETECTED.
!           *** EXTEND LAST DEPTH BEYOND MAX RANGE.
            TRACK(NTRK+1,1,L)=1.0E+38
            TRACK(NTRK+1,2,L)=TRACK(NTRK,2,L)
            R2(L)=TRACK(1,1,L)
            Z2(L)=TRACK(1,2,L)
          ENDIF
   24   CONTINUE
!
      ELSE
!       *** USER WRITTEN SUBROUTINE UBOTTOM SUPPLIES PROFILE.
        CALL UBOTTOM
      ENDIF
!
!     *** FIND BOTTOM SEGMENT WHICH CONTAINS STARTING RANGE
      DO 31 L=1,NSOL
        ITRK(L)=1
   25      R1(L)=R2(L)
        Z1(L)=Z2(L)
        ITRK(L)=ITRK(L)+1
!
        IF (ITRK(L) > MXTRK) THEN
          ITEMP=MXTRK
          WRITE(NPU,29) ITEMP
   29     FORMAT(1X,'ERROR. BOTTOM PROFILE MISSING OR ARRAY OVERFLOW. ',    &
          'MAX IS ',I5, ',. CHECK PARAMETER MXTRK.')
          STOP
        ELSE
        ENDIF
!
        R2(L)=TRACK(ITRK(L),1,L)
        Z2(L)=TRACK(ITRK(L),2,L)
!
!       *** ADVANCE TRACK IF NECESSARY SO THAT R1 <= RA < R2
        IF (RA >= R2(L)) GO TO 25
!
!       *** BE SURE THAT R1 <= RA
        IF (R1(L) > RA) THEN
          WRITE(NPU,*)'DEPTH OF BOTTOM AT INITIAL RANGE MISSING,'
          WRITE(NPU,*)'OR PARAMETER MXTRK SHOULD BE INCREASED'
          STOP
        ELSE
        ENDIF
!
!       *** COMPUTE SLOPE OF BOTTOM IN RANGE
        THETA(L)=ATAN2(Z2(L)-Z1(L),R2(L)-R1(L))
!
!       *** END OF DO LOOP
   31 CONTINUE
!
!     *** READ RANGE OF SVP
      READ(NIU,*,END=33) RSVP
!
!     *** SVP BEYOND START RANGE?
      IF (RSVP > RA) GO TO 33
!
!     *** NO. DETERMINE IF SVP IN RUNSTREAM OR SUPPLIED BY SUBROUTINE USVP3D.
   32 READ(NIU,*,END=33) KSVP
      IF (KSVP == 0) CALL SVP3D
      IF (KSVP /= 0) CALL USVP3D
!
!     *** ERROR IN SVP?
      IF (NSVP /= 0) GO TO 35
!     *** YES.
   33 WRITE(NPU,34) RA
   34 FORMAT(1X,'SVP MISSING OR INPUT ERROR. RANGE = ',F8.1,' M.')
      STOP
!
   35 CONTINUE
!     *** NO ERROR DETECTED.
!     *** IF C0 NOT SPECIFIED, SET C0 TO AVERAGE SPEED OF FIRST PROFILES
!     *** IN FIRST LAYERS AT INITIAL RANGE
!
      IF (C0 == 0.0) THEN
        ZWSVP=0.0
        DO L=1,NSOL
          NWSVP=IXSVP(1,L)
          ZWSVP=ZWSVP+ZSVP(NWSVP,L)
          DO I=2,NWSVP
            C0=C0+(ZSVP(I,L)-ZSVP(I-1,L))*(CSVP(I-1,L)+.5*(CSVP(I,L)-    &
            CSVP(I-1,L)))
          ENDDO
        ENDDO
        C0=C0/ZWSVP
      ELSE
      ENDIF
!
!     *** COMPUTE REFERENCE WAVE NUMBER
      XK0=2.0*PI*FRQ/C0
!
!     *** COMPUTE ATTENUATION - SACLANT MEMO SM-121 (JENSEN + FERLA)
!     *** MODIFIED AS FOLLOWS:
!     *** IF INPUTTED BETA IS LT 0.0, ALPHA IS COMPUTED IN DB/METER
!     *** AND USED TO COMPUTE BETA IN DB/WAVELENGTH ROUTINE INDX3D.
!
      ALPHA=FRQ*FRQ*(.007+(.155*1.7)/(1.7*1.7+FRQ*FRQ*.000001))*1.0E-09
!
!     *** ADJUST LAYER DEPTHS IN CASE BOTTOM SLOPES AND RA > RSVP.
!     *** ASSUMES LAYERS FOLLOW BOTTOM CONTOUR.
!
      DO L=1,NSOL
        NLYR=NLYRS(L)
        DO ILYR=1,NLYR
          ZLYR(ILYR,L)=ZLYR(ILYR,L)+(RA-RSVP)*TAN(THETA(L))
        ENDDO
      ENDDO
!
!     *** GET RANGE OF NEXT SVP
!     *** IF END, NO MORE PROFILES FOLLOW.
!
   53 READ(NIU,*,END=55) RSVP
!
!     *** TEST FOR 3D INPUT RUNSTREAM - USED DURING DEBUG PHASE.
!     *** WHEN COMPARING 2D SOLUTION WITH AZIMUTH INDEPENDENT NX2D SOLUTIONS.
!     *** SEE NOTE 2 REGARDING -1,-1 IN DESCRIPTION OF INPUT PARAMETERS.
!
      IF (NDIM == 1 .AND. RSVP < 0.0)GO TO 53
      GO TO 56
!
!     *** ONLY ONE SVP - SET RSVP LARGE SO SAME SVP USED FOR ENTIRE PROBLEM
   55 RSVP=1.0E+38
!
   56 CONTINUE
!
!     *** IF STARTING FIELD IS BEYOND NEXT SVP, GO BACK AND GET NEXT SVP
      IF (RSVP <= RA) GO TO 32
!
!     *** ABSORBING LAYER REQUESTED?
      IF (ITYPEB == 3) THEN
!       *** YES.
!       *** EXTEND BOTTOM 4/3 MAX DEPTH IF ZA = ZERO
        IF (ZA == 0.0) ZA=4.0*ZBOT/3.0
!       *** ABSORBING LAYER DEEP ENOUGH?
        IF (ZA < ZBOT) THEN
!         *** NO. EXTEND BOTTOM 4/3 MAX DEPTH.
          WRITE(NPU,57)
   57     FORMAT(1X,'ERROR. ZA RESET TO MAX DEPTH OF BOTTOM LAYER.')
          ZA=4.0*ZBOT/3.0
        ELSE
        ENDIF
!
!       *** INSERT PARAMETERS FOR EXTENDED BOTTOM IN APPROPRIATE ARRAYS
        ZMAX=ZA
        DO L=1,NSOL
          NLYR=NLYRS(L)
          ZG=ZLYR(NLYR,L)
          NSVP=IXSVP(NLYR,L)
          IF (NLYR > 1)ZG=ZG-ZLYR(NLYR-1,L)
          NLYRS(L)=NLYRS(L)+1
          NLYR=NLYRS(L)
!         *** store depth of artificial absorbing layer
          ZLYR(NLYR,L)=ZA
!         *** use bottom density and attenuation of above layer
          RHO(NLYR,L)=RHO(NLYR-1,L)+ZG*RHOG(NLYR-1,L)
          RHOG(NLYR,L)=0.0
          BETA(NLYR,L)=BETA(NLYR-1,L)+ZG*BETAG(NLYR-1,L)
          BETAG(NLYR,L)=0.0
!         *** use bottom speed of above layer
          IXSVP(NLYR,L)=NSVP+2
          ZSVP(NSVP+1,L)=ZSVP(NSVP,L)
          CSVP(NSVP+1,L)=CSVP(NSVP,L)
          ZSVP(NSVP+2,L)=ZLYR(NLYR,L)
          CSVP(NSVP+2,L)=CSVP(NSVP+1,L)
        ENDDO
      ELSE
!       *** IF BOTTOM NOT EXTENDED AND ZA=0, SET ZA TO MAX DEPTH INPUTTED
!       *** WITH FIRST SVP. ASSUMES DEPTH IS CONSTANT IN AZIMUTH.
        IF (ZA == 0.0) ZA=ZLYR(NLYR,1)
      ENDIF
!
!     *** IF N NOT SPECIFIED, SET N SO THAT DZ  <=  1/10 WAVELENGTH
      IF (N == 0) THEN
        N=INT((10*ZA*FRQ/C0)+1)
      ELSE
      ENDIF
!
!     *** COMPUTE RECEIVER DEPTH INCREMENT. DZ MAY BE SUCH THAT RECEIVERS DO
!     *** NOT LIE EXACTLY ON LAYER INTERFACES
      DZ=ZA/N
!
!     ADJUST N SO THAT DEPTH (N+1)*DZ IS BOTTOM
      IF (ITYPEB == 1 .OR. ITYPEB == 3)N=N-1
!
!     *** COMPUTE TOTAL SOLUTION DEPTHS
      NTOT=N*NSOL
!
!     *** IF RANGE STEP NOT SPECIFIED, SET IT EQUAL TO 1/10 WAVELENGTH.
!     *** THIS STEP SIZE IS ARBITRARY.
!
      IF (DR == 0.0) DR=.1*C0/FRQ

!WYX BEG ------------------------------------------------------
!
      ALLOCATE( D(N, NSOL) )
      ALLOCATE( U(0:N+1, 0:NSOL+1) )

      ALLOCATE( AL(N, NSOL) )
      ALLOCATE( AU(N, NSOL) )
      ALLOCATE( AM(N, NSOL) )
      ALLOCATE( BL(N, NSOL) )
      ALLOCATE( BU(N, NSOL) )
      ALLOCATE( BM(N, NSOL) )

      ALLOCATE( XN1(N, NSOL) )
      ALLOCATE( XN2(N, NSOL) )
      ALLOCATE( RHO1(N, NSOL) )
      ALLOCATE( RHO2(N, NSOL) )


      ALLOCATE( DP(1:N) )
      ALLOCATE( DS(1:N) )
      ALLOCATE( PORTX(1:N) )
      ALLOCATE( PORTY(1:N) )
      ALLOCATE( STBDX(1:N) )
      ALLOCATE( STBDY(1:N) )

      ALLOCATE( BOTX(1:NSOL+2) )
      ALLOCATE( BOTY(1:NSOL+2) )
      ALLOCATE( SURX(1:NSOL+2) )
      ALLOCATE( SURY(1:NSOL+2) )

!
!WYX END ------------------------------------------------------

!
!     *** GET STARTING FIELD
      CALL SFLD3D
!
!     *** COMPUTE DEPTH WRITE INCREMENT TO NEAREST DZ
      IF (WZ1 < DZ)WZ1=DZ
      IWZ1=INT(WZ1/DZ+.5)
!
      IF (WZ2 < DZ)WZ2=DZ
      IF (ITYPEB == 1 .OR. ITYPEB == 3 .AND. WZ2 > N*DZ)WZ2=N*DZ
       IWZ2=INT(WZ2/DZ+.5)
!
      IF (WDZ < DZ)WDZ=WZ1
      IWZ=INT(WDZ/DZ+.5)
      WDZ=IWZ*DZ
!
      NZ=(IWZ2-IWZ1)/IWZ+1
!
!     *** COMPUTE DEPTH PRINT INCREMENT TO NEAREST DZ
      IPZ=INT(PDZ/DZ+.5)
      IF (PDZ > 0.0 .AND. IPZ == 0) IPZ=1
      PDZ=IPZ*DZ
!
!     *** PRINT INPUT PARAMETERS
#ifdef NPU_ENABLE
      CALL PRINTP
#endif
!
   84 FORMAT(3X,'LAYER ',I3,3X,'DEPTH = ',F8.1,3X,    &
      'DENSITY = ',F6.2,3X,'GRADIENT = ',F7.4,/,3X,'BETA = ',F8.3,    &
      3X,'GRADIENT = ',F7.4)
   93 FORMAT(//1X,'SOUND VELOCITY PROFILE AT RANGE ',F10.1,' METERS',/)
   94 FORMAT(1X,I4,3X,F8.1,3X,F8.1)
#ifdef NOU_ENABLE
!
!     *** OUTPUT REQUESTED?
      IF (WDR /= 0.0) THEN
!       *** YES. WRITE SELECTED PARAMETERS TO NPU_ENABLE FILE.
!       CALL ASSIGN(NOU,'FOR3D.OUT')
!ljh
        OPEN(UNIT=NOU,FILE='FOR3D.OUT')
        WRITE(NOU,994)NDIM,FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,    &
        ITYPPW,ITYPSW,FLDW,NSEC,NSOL,RMAX,DR,WDR,WZ1,WZ2,WDZ,DZ,    &
        DOUGRA,NDIV,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
  994   FORMAT(1X,I1,3F10.2,I3,2F8.0,I5,5I2,F9.3,/,2I5,F10.1,F5.1,    &
        5F7.1,F10.1,I4,/,(5F15.7))
!
!       *** WRITE STARTING FIELD TO OUTPUT FILE.
        DO 97 J=1,NSOL,NSOLW
          IF (NDIM /= 1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC)
          WRITE(NOU,995) ANG,NZ,RA,WZ1,WDZ,(U(I,J),I=IWZ1,IWZ2,IWZ)
  995     FORMAT(1X,F10.5,I6,F8.0,F8.1,F8.1,/,2(E15.7,E15.7))
   97   CONTINUE
      ELSE
!       *** NO NPU_ENABLE REQUESTED.
      ENDIF
#endif
!
!
!cljh
!
!     *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE PRINTED
      XPR=RA+PDR
!
!     *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE RECORDED
      XWR=RA+WDR
      IF (WDR == 0.0) XWR=RA+RMAX+1.0
!
!     *** SAVE RANGE STEP
      OLDR=DR
!
      N1=1
!
!     ******************************************************************
!     ******************************************************************
!     *** MAIN LOOP STARTS HERE                                      ***
!     *** SOLUTION WILL BE ADVANCED FROM RANGE RA TO RANGE RA+DR     ***
!     ******************************************************************
!     ******************************************************************
!
      MAIN_LOOP: DO WHILE (.TRUE.)
      DRLAST=DR
!
!     *** Compute index of refraction table at present range
      CALL INDX3D
!
!     *** This section of code determines whether or not to restore step
!     *** size to the original DR.
!
!ljh  Well, this looks important!  How far does "Next Section" extend???
!     *** Comment next section if input DR is to be used
      IF (ALL(THETA(1:NSOL) == 0.0)) THEN
!       *** flat bottom - Restore DR if it has been changed
        DR=OLDR
      ENDIF

      DO L=1,NSOL
        DRMIN(L)=DR
!       *** Does SVP change before bottom profile?
        IF (RSVP < R2(L)) THEN
!         *** Yes. Does SVP change before next solution range?
          IF (RA+DRMIN(L) > RSVP) THEN
!           *** Yes. Adjust DR so that solution advances to RSVP
            DRMIN(L)=RSVP-RA
          ENDIF
          CYCLE
        ENDIF
!       *** bottom profile changes before or at same range as SVP
!
!       *** DETERMINE IF BOTTOM DEPTHS ARE TO BE UPDATED
!       *** FIRST PASS - RA WILL BE  <  R2
!
!       *** UPDATE BOTTOM DEPTHS?
        IF (RA < R2(L)) THEN
!         *** NO.
!         *** DETERMINE IF RANGE STEP TOO LARGE .
          IF (RA+DRMIN(L) > R2(L)) THEN
!           *** RANGE STEP IS TOO LARGE - RESET DR - ADV SOLUTION TO R2
            DRMIN(L)=R2(L)-RA
          ENDIF
          CYCLE
        ENDIF
!
!       *** UPDATE BOTTOM DEPTHS
        DO WHILE(.TRUE.)
          R1(L)=R2(L)
          Z1(L)=Z2(L)
          ITRK(L)=ITRK(L)+1
          R2(L)=TRACK(ITRK(L),1,L)
          Z2(L)=TRACK(ITRK(L),2,L)
!
!         *** TWO DEPTHS AT SAME RANGE INDICATE VERTICAL DISCONTINUITY.
!         *** ADVANCE TRACK FORWARD.
          IF (R1(L) /= R2(L)) EXIT
        ENDDO
!
!       *** Restore DR in case bottom is flat
        DRMIN(L)=OLDR
!
!       *** Compute slope of bottom
        THETA(L)=ATAN2(Z2(L)-Z1(L),R2(L)-R1(L))
!
!       *** IF BOTTOM IS NOT FLAT, COMPUTE NEW RANGE STEP
        IF (THETA(L) /= 0) THEN
!         *** COMPUTE NEW RANGE STEP. DR MAY BECOME VERY SMALL.
          DRMIN(L)=ABS(DZ*COS(THETA(L))/SIN(THETA(L)))
        ENDIF
!       *** IF RANGE STEP TOO LARGE, PRINT WARNING MESSAGE. RECOMPUTE DR.
        IF (RA+DRMIN(L) > R2(L)) THEN
!         *** ADJUST DR AND CONTINUE.
          DRMIN(L)=R2(L)-RA
        ENDIF
        IF (RA+DRMIN(L) > RSVP) DRMIN(L)=RSVP-RA
!
!       *** PRINT BOTTOM DEPTHS
#ifdef NPU_ENABLE
        IF (IBOT /= 0) THEN
          WRITE(NPU,93) RA
          IF (NDIM /= 1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC)
          WRITE(NPU,*)' RELATIVE BEARING =  ',ANG,' DEGREES'
          TH=THETA(L)*DEG
          WRITE(NPU,123) R1(L),Z1(L),R2(L),Z2(L),TH
  123     FORMAT(//1X,'BOTTOM PROFILE ',//,1X,    &
          'R1 = ',F10.1,' M',/,1X,    &
          'Z1 = ',F8.1,' M',/,1X,    &
          'R2 = ',F10.1,' M',/,1X,    &
          'Z2 = ',F8.1,' M',/,1X,    &
          'SLOPE = ',F8.1,' DEG')
        ENDIF
#endif
      ENDDO
!     *** END OF LOOP
!
!     *** Select smallest DR within following arbitrary limits.
!     *** Maximum allowed is input DR.
!     *** Minimum allowed is 0.5*(DZ*DZ).
      DR = MINVAL(DRMIN(1:NSOL))
      DR = MAX(DR, 0.5*DZ*DZ)
      !DR = MIN(DR, OLDR)
      IF (DR > OLDR) DR = OLDR

      IF (DR /= DRLAST)WRITE(NPU,*)'AT RANGE ',RA,' DR WAS CHANGED FROM ',DRLAST,' TO ',DR
!
!     *** Compute matrix A at present range
      CALL AMIFD3
!
!     *** If 3D run, compute B matrix.
      IF (NDIM == 3) CALL BMIFD3
!
!     *** Advance range one step
      RA=RA+DR
      NOLD=N
!
!     *** Compute right hand side
      CALL RHS
!
!     *** March solution 1 DR forward
      CALL TWOSTEP
!
!     *** APPLY ABSORPTION IF ITYPEB = 3
      IF (ITYPEB == 3) THEN
        DO 180 J=1,NSOL
          NLYR=NLYRS(J)
          MM=INT(ZLYR(NLYR-1,J)/DZ)
          NA=N-MM
          IF (NA > 0) GO TO 175
          IF (NA >= -1) GO TO 180
          WRITE(NPU,174) RA
  174     FORMAT(1X,'ERR IN THICKNESS OF ABSORBING LAYER AT ',F8.1,' M')
          STOP
!
  175     CONTINUE
          DO 179 I=1,NA
            ATT=EXP(-.01*DR*EXP(-((I-NA)/(NA/3.0))**2))
            U(MM+I,J)=U(MM+I,J)*ATT
            IF (J > 1)GO TO 176
            IF (NDIM == 3 .AND. ITYPPW == 2)PORTY(MM+I)=PORTY(MM+I)*ATT
  176       IF (J < NSOL)GO TO 179
            IF (NDIM == 3 .AND. ITYPSW == 2)STBDY(MM+I)=STBDY(MM+I)*ATT
  179     CONTINUE
  180   CONTINUE
      ELSE
      ENDIF
!
#ifdef NOU_ENABLE
!     *** IF SOLUTION IS TO BE WRITTEN ON DISK,
!     *** WRITE SELECTED PRESSURE FIELD AT RANGE RA
      IF (XWR <= RA) THEN
        IF (WZ1 < DZ) WZ1=DZ
        IWZ1=NINT(WZ1/DZ)
!
        IF (WZ2 < DZ) WZ2=DZ
        IF (ITYPEB == 1 .OR. ITYPEB == 3 .AND. WZ2 > N*DZ) WZ2=N*DZ
        IWZ2=NINT(WZ2/DZ)
        IF (WDZ < DZ) WDZ=WZ1
!
        IWZ=NINT(WDZ/DZ)
        WDZ=IWZ*DZ
!
        NZ=(IWZ2-IWZ1)/IWZ+1

        DO J=1,NSOL,NSOLW
          IF (NDIM /= 1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC)
          WRITE(NOU,995) ANG,NZ,RA,WZ1,WDZ,(U(I,J),I=IWZ1,IWZ2,IWZ)
        ENDDO
!       *** DETERMINE NEXT RANGE AT WHICH TO WRITE SOLUTION ON DISK
        XWR=XWR+WDR
      ENDIF
#endif
!
!     *** DETERMINE IF SOLUTION IS TO BE PRINTED
      IF (XPR <= RA .AND. IPZ /= 0 .AND. PDR /= 0.0) THEN
!         *** PRINT RANGE
        RTEMP=RA/1000.0
!
#ifdef NOU_ENABLE
        OPEN(UNIT=NOU,FILE='FOR3D.OUT')
#endif
#ifdef NPU_ENABLE
!       *** Compute and print propagation loss at each IPZ'th depth
  275   FORMAT(2X,'RECEIVER',4X,'DEPTH(M)',4X,'LOSS(DB)',14X,'U(I)')
        JSOL=NSOL/2
        DO J=1,NSOL,NSOLP
          IF (NDIM /= 1)ANG=-FLDW*PI*(.5*NSEC-(J-1))/(180.0*NSEC)
          WRITE(NPU,277)J,RA,ANG*DEG
  277     FORMAT(/,2X,'SOLUTION',I4,',  RANGE = ',F8.1,' M',    &
          ',  RELATIVE BEARING = ',F6.1,    &
          ' DEGREES',/)
          WRITE(NPU,275)
          UA=0
          DO I=IPZ,N,IPZ
            ZI=I*DZ
            UA=U(I,J)
            IF (IHNK == 1)UA=UA*HNKL(XK0*RA)
            PL=ABS(UA)
            IF (PL <= 0.0) CYCLE
            PL=-20.0*DLOG10(PL)                    ! dlog10(real*8)
            IF (IHNK == 0) PL=PL+10.0*ALOG10(RA)
            WRITE(NPU,294) I,ZI,PL,UA
  294       FORMAT(2X,I5,(3X,F10.2,3X,F10.3),3X,'(',E12.5,2X,E12.5,' )')
!ljh        Another important bit of code to change, hidden deep in the program!
!           *** UNCOMMENT IF EXACT SOLUTION SUBROUTINE AVAILABLE.
!           CALL UEXACT
!           WRITE(NPU,*)' '
          ENDDO
        ENDDO
#endif
!
!       *** DETERMINE NEXT RANGE AT WHICH TO PRINT SOLUTION
        XPR=XPR+PDR
      ENDIF
!
!     *** SWITCH TO DOUGLAS METHOD?
      IF (DOUGRA /= 0.0 .AND. RA >= DOUGRA) THEN
!       *** YES.
        OLDZ=DZ
        NOLD=N
        NTOLD=NTOT
        IF (ITYPEB == 1) N=N+1
        IF (ITYPEB == 3) N=N+1
!
!       *** DIVIDE N BY NDIV AND MAKE NECESSARY ADJUSTMENTS.
        N=INT(N/NDIV+0.5)
        DZ=ZA/N
        IF (ITYPEB == 1)N=N-1
        IF (ITYPEB == 3)N=N-1
        NTOT=N*NSOL
        IPZ=NINT(PDZ/DZ)
!
!       *** ADJUST SOLUTION FIELD.
IF (NOLD /= N) WRITE(*,*) "ERROR ADDED BY WYX: NOLD /= N. ", NOLD, N
        DO J=1,NSOL
          MOLD=(J-1)*NOLD
          M=(J-1)*N
          DO I=1,N
            ZI=I*DZ
!           *** AVOID ROUNDOFF PROBLEM.
            K=NINT(ZI/OLDZ)
            U(I,J)=U(K,J)+(ZI-K*OLDZ)*(U(K+1,J)-U(K,J))/OLDZ
          ENDDO
        ENDDO
!
!       *** ADJUST WALL CONDITIONS.
        DO I=1,N
          ZI=I*DZ
!         AVOID ROUNDOFF PROBLEM.
          K=INT(ZI/OLDZ+.5)
          PORTY(I)=PORTY(K)
          STBDY(I)=STBDY(K)
        ENDDO
!
!       *** SET ALFA FOR DOUGLAS METHOD.
        ALFA=1.0/12.0
!       *** SET DOUGRA LARGE SO THAT WE DON'T ENTER THIS CODE AGAIN.
        DOUGRA=1.0E+38
!
#ifdef NPU_ENABLE
!       WRITE(NPU,*)' '
        WRITE(NPU,*)'AT RANGE ',RA,' METERS. SWITCH FROM CRANK-NICOLSON TO DOUGLAS METHOD.'
        WRITE(NPU,*)'CRANK-NICOLSON    N ,DZ, NTOT ',NOLD,OLDZ,NTOLD
        WRITE(NPU,*)'DOUGLAS           N ,DZ, NTOT ',N,     DZ,NTOT
#endif
!CCC        CALL SECOND(CPT2)
!ljh        WRITE(NPU,370)CPT2-CPT1
!ljh  370   FORMAT(1X,'ELAPSED TIME = ',F12.2,' CRAY XMP 2/8 CPU SECONDS.')
!       WRITE(NPU,*)' '
!
!       *** RECOMPUTE DEPTH WRITE INCREMENT TO NEAREST DZ
        IF (WZ1 < DZ)WZ1=DZ
        IWZ1=NINT(WZ1/DZ)
!
        IF (WZ2 < DZ)WZ2=DZ
        IWZ2=NINT(WZ2/DZ)
!
        IF (WDZ < DZ)WDZ=WZ1
        IWZ=NINT(WDZ/DZ)
        WDZ=IWZ*DZ
!
        NZ=(IWZ2-IWZ1)/IWZ+1
!
!       *** RECOMPUTE DEPTH PRINT INCREMENT TO NEAREST DZ
        IPZ=NINT(PDZ/DZ)
        IF (PDZ > 0.0 .AND. IPZ == 0) IPZ=1
        PDZ=IPZ*DZ
!
      ELSE
      ENDIF      ! End of If SWITCH TO DOUGLAS METHOD?, line 1211
!
!     *** Terminate run if solution at maximum range has been obtained
      IF (RA >= RMAX) EXIT MAIN_LOOP
!
!     *** Get environment for next range
!     *** Read new SVP profile flag
      IF (RA >= RSVP) READ(NIU,*,END=33) KSVP
!
!     *** KSVP MAY HAVE BEEN SET IN USER ROUTINE USVP3D.
!     *** IF KSVP NOT 0, SUBROUTINE USVP3D IS CALLED. USER IS
!     *** RESPONSIBLE FOR ADJUSTING DEPTHS OF LAYERS.
      IF (KSVP /= 0) THEN
        CALL USVP3D
      ELSE
!       *** NEW SVP AT ADVANCED RANGE?
        IF (RA < RSVP) THEN
!         *** UPDATE DEPTHS OF LAYERS
!         *** ASSUMES LAYERS ARE PARALLEL AND FOLLOW BOTTOM CONTOUR
!         *** IF ITYPEB = 3, BOTTOM OF ARTIFICIAL LAYER REMAINS FLAT.
          DO L=1,NSOL
            DZZ=DR*TAN(THETA(L))
!           *** IS BOTTOM FLAT?
            IF (THETA(L) == 0.0) CYCLE
            NLYR=NLYRS(L)
            MLYR=NLYRS(L)
            IF (ITYPEB == 3) MLYR=NLYRS(L)-1
            DO ILYR=1,MLYR
              ZLYR(ILYR,L)=ZLYR(ILYR,L)+DZZ
              IF (ITYPEB == 3 .AND. ZLYR(NLYR,L) < ZLYR(MLYR,L)) THEN
                WRITE(NPU,*)'ERROR. ARTIFICIAL LAYER NOT DEEP ENOUGH.'
                STOP
              ELSE
              ENDIF
            ENDDO
          ENDDO
!
!         *** ADJUST DEPTHS OF PROFILES IN SLOPING LAYERS
!         *** ASSUMES SAME SVP IN LAYERS
          DO L=1,NSOL
            DZZ=DR*TAN(THETA(L))
            NLYR=NLYRS(L)
            NWSVP=IXSVP(1,L)
            NSVP=IXSVP(NLYR,L)
            DO I=NWSVP+1,NSVP
              ZSVP(I,L)=ZSVP(I,L)+DZZ
            ENDDO
          ENDDO
          GO TO 570
        ENDIF

!       *** GET NEXT SVP FROM INPUT RUNSTREAM.
        CALL SVP3D
      ENDIF

!     *** ERROR DETECTED?
      IF (NSVP == 0) GO TO 33
!
!     *** NO. READ RANGE OF NEXT PROFILE?
      IF (RA >= RSVP) THEN
        READ(NIU,*,END=549) RSVP
!       *** CHECK FOR 3D RUNSTREAM
        IF (NDIM == 1 .AND. RSVP < 0.0)READ(NIU,*,END=549) RSVP
      ELSE
      ENDIF
      GO TO 550
!
  549 CONTINUE
!     *** SET RSVP LARGE SO THAT LAST PROFILE IS USED FOR REMAINDER OF PROBLEM
      RSVP=1.0E+38
!
  550 CONTINUE
!     *** ARTIFICIAL ABSORBING LAYER?
      IF (ITYPEB == 3) THEN
!       *** YES. UPDATE DENSITY, ATTEN AND SPEED.
        DO L=1,NSOL
          NLYR=NLYRS(L)
          ZG=ZLYR(NLYR,L)
          NSVP=IXSVP(NLYR,L)
          IF (ZA < ZLYR(NLYR,L)) CYCLE
          IF (NLYR > 1)ZG=ZLYR(NLYR,L)-ZLYR(NLYR-1,L)
          NLYRS(L)=NLYRS(L)+1
          NLYR=NLYRS(L)
          ZLYR(NLYR,L)=ZA
          RHO(NLYR,L)=RHO(NLYR-1,L)+ZG*RHOG(NLYR-1,L)
          RHOG(NLYR,L)=0.0
          BETA(NLYR,L)=BETA(NLYR-1,L)+ZG*BETAG(NLYR-1,L)
          BETAG(NLYR,L)=0.0
          IXSVP(NLYR,L)=NSVP+2
          ZSVP(NSVP+1,L)=ZSVP(NSVP,L)
          CSVP(NSVP+1,L)=CSVP(NSVP,L)
          ZSVP(NSVP+2,L)=ZLYR(NLYR,L)
          CSVP(NSVP+2,L)=CSVP(NSVP+1,L)
        ENDDO
      ELSE
      ENDIF
#ifdef NPU_ENABLE
!
!     *** PRINT SVP
      IF (ISVP /= 0) THEN
        WRITE(NPU,93) RA
        DO L=1,NSOL
          IF (NDIM /= 1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC)
          WRITE(NPU,*)' RELATIVE BEARING =  ',ANG,' DEGREES'
          K=0
          NLYR=NLYRS(L)
          DO ILYR=1,NLYR
            WRITE(NPU,84)ILYR,ZLYR(ILYR,L),RHO(ILYR,L),RHOG(ILYR,L),    &
            BETA(ILYR,L),BETAG(ILYR,L)
            J=K+1
            K=IXSVP(ILYR,L)
            DO I=J,K
              WRITE(NPU,94) I,ZSVP(I,L),CSVP(I,L)
            ENDDO
          ENDDO
        ENDDO
      ELSE
      ENDIF
#endif
  570 CONTINUE

      ENDDO MAIN_LOOP


!     *** TERMINATE RUN
#ifdef NOU_ENABLE
      IF (WDR /= 0.0) CLOSE(NOU)
#endif
END PROGRAM
!                                                             **********
!                                                             * AMIFD3 *
!                                                             **********
SUBROUTINE AMIFD3
!     ******************************************************************
!     ******************************************************************
!     * COMPUTE DIAGONALS FOR 'A' MATRIX                               *
!     ******************************************************************
      USE commons, ONLY: AU, AL, AM, DR, DZ, N, NSOL, XK0, ALFA, XN1, XN2, &
                         RHO1, RHO2
      USE omp_lib
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, L, M
      REAL    (KIND=8) :: R12, RSUM, T2
      COMPLEX          :: T1, T3, XN12, XNSUM

!     *** RA has been incremented by DR in main program
      T1=CMPLX(.25,-.25*XK0*DR)
      T2=1.0/(XK0*XK0*DZ*DZ)
      T3=T1*T2
      
       !$OMP PARALLEL DO PRIVATE(M,XN12,R12,RSUM,XNSUM)
      DO L=1,NSOL
        M=(L-1)*N
        DO I=1,N
          XN12=T1*((XN1(I,L)-1)/RHO1(I,L)+(XN2(I,L)-1)/RHO2(I,L))
          R12=1.0/RHO1(I,L)+1.0/RHO2(I,L)
          RSUM=RHO1(I,L)+RHO2(I,L)
          XNSUM=T1*(XN1(I,L)+XN2(I,L)-2.0)
!**       UPPER DIAGONAL
          AU(I,L)=(T3+ALFA+ALFA*T1*(XN2(I,L)-1))/RHO2(I,L)
          IF (I == N) AU(I,L)=0.0
!**       MAIN DIAGONAL
          AM(I,L)=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!**       LOWER DIAGONAL
          AL(I,L)=(T3+ALFA+ALFA*T1*(XN1(I,L)-1))/RHO1(I,L)
          IF (I == 1)AL(I,L)=0.0
        ENDDO
      ENDDO
END SUBROUTINE
!                                                             **********
!                                                             * BCON3D *
!                                                             **********
SUBROUTINE BCON3D
!     *****************************************************************
!     *** BOTTOM BOUNDARY CONDITION SUBROUTINE
!     *****************************************************************
      USE commons, ONLY: N, ITYPEB, NSOL, NPU, BOTX, BOTY, U
      IMPLICIT NONE
      INTEGER (KIND=4) :: IERR, L

      SELECT CASE (ITYPEB)

!     *** PRESSURE RELEASE BOTTOM
      CASE (0)
        DO L=1,NSOL+2
          BOTY(L)=0.0
          BOTX(L)=0.0
        ENDDO

!     *** USER BOTTOM BOUNDARY CONDITION - USER WRITES SUBROUTINE UBCON3D
      CASE (1)
        CALL UBCON3D

!     *** ARTIFICIAL ABSORBING BOTTOM - FLAT
      CASE (3)
        DO L=1,NSOL+2
          BOTY(L)=0.0
          BOTX(L)=0.0
        ENDDO

      CASE DEFAULT
        WRITE(NPU,50)ITYPEB
   50   FORMAT(1X,'ERROR - BOTTOM BOUNDARY CONDITION ',I2,    &
        ' NOT IMPLEMENTED')
        IERR=1

      END SELECT

      U(N+1, 0:NSOL+1) = BOTY(1:NSOL+2)
END SUBROUTINE
!                                                             **********
!                                                             * BMIFD3 *
!                                                             **********
SUBROUTINE BMIFD3
!     ******************************************************************
!     * COMPUTE DIAGONALS FOR B MATRIX                                 *
!     ******************************************************************
      USE commons, ONLY: DR, N, NSOL, PHI, RA, XK0, RHO1, &
                         BL, BU, BM
      USE omp_lib
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, J, L, LM1, LP1, M
      REAL    (KIND=8) :: R34, RH3, RH4
      COMPLEX          :: T1

      T1=-CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA+.5*DR)*(RA+.5*DR)*PHI*PHI))
      
       !$OMP PARALLEL DO PRIVATE(M,LM1,LP1,J,RH3,RH4,R34)
      DO I=1,N
        M=(I-1)*NSOL
        DO L=1,NSOL
          LM1=L-1
          IF (LM1 == 0)LM1=1
          LP1=L
          IF (LP1 > NSOL)LP1=NSOL
          J=(LM1-1)*N
          RH3=RHO1(I,LM1)
          J=(LP1-1)*N
          RH4=RHO1(I,LP1)
          R34=(1.0/RH3+1.0/RH4)
!         *** UPPER DIAGONAL
          BU(I,L)=.5*(rh3+rh4)*(T1)/RH4
          IF (L == NSOL)BU(I,L)=0.0
!         *** MAIN DIAGONAL
          BM(I,L)=1.0+(-T1)*.5*(rh3+rh4)*R34
!         *** LOWER DIAGONAL
          BL(I,L)=.5*(rh3+rh4)*(T1)/RH3
          IF (L == 1)BL(I,L)=0.0
        ENDDO
      ENDDO
END SUBROUTINE
!                                                             **********
!                                                             *  HNKL  *
!                                                             **********
COMPLEX FUNCTION HNKL(X)
!     ******************************************************************
!     * HANKEL FUNCTION H0(1) - POLYNOMIAL APPROXIMATION               *
!     * HANDBOOK OF MATH FUNCTIONS - N.B.S. - NOV 1967                 *
!     ******************************************************************
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, PARAMETER :: PI = 3.14159265358979
      REAL :: J0, F0, T0, Y, Y0
!
      IF (X > 3.) THEN
!
!     *** (3.0 <= X < INFINITY)
      Y=3.0/X
      F0=  0.79788456+Y*(-0.00000077+Y*(-0.00552740+Y*(-0.00009512    &
                     +Y*(+0.00137237+Y*(-0.00072805+Y*(+0.00014476))))))
      T0=X-0.78539816+Y*(-0.04166397+Y*(-0.00003954+Y*(+0.00262573    &
                     +Y*(-0.00054125+Y*(-0.00029333+Y*(+0.00013558))))))
      HNKL=F0*CEXP(CMPLX(0.0,T0))/SQRT(X)

      ELSE
!
!     *** (-3.0 <= X <= 3.0)
      Y=X*X/9.0
      J0=1.+Y*(-2.2499997+Y*(+1.2656208+Y*(-0.3163866+Y*(+0.0444479    &
           +Y*(-0.0039444+Y*(+0.0002100))))))
!
!     *** (0.0 < X <= 3.0)
      Y0=2.0*LOG(0.5*X)*J0/PI+0.36746691    &
           +Y*(+0.60559366+Y*(-0.74350384+Y*(+0.25300117+Y*(-0.04261214    &
           +Y*(+0.00427916+Y*(-0.00024846))))))
      HNKL=CMPLX(J0,Y0)

      ENDIF
END FUNCTION
!                                                             **********
!                                                             * INDX3D *
!                                                             **********
SUBROUTINE INDX3D
!     ******************************************************************
!     * Compute index of refraction table
!     ******************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons, ONLY: ALPHA, C0, DZ, FRQ, N, NLYR, NPU, NSOL, ZI, &
                         NLYRS, IXSVP, ZLYR, IXSVP, ZSVP, CSVP, RHO, &
                         RHOG, BETA, BETAG, &
                         RHO1, RHO2, XN1, XN2
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, IEX, ILYR, ISOL, J, L, M, MM
      REAL    (KIND=8) :: BETA1, BETA1G, BETA2, BETA2G, C1, C2, CI, &
                          RH1, RH1G, RH2, RH2G, XN, ZG1, ZG2

      IEX = 0
      DO ISOL=1,NSOL
        MM=(ISOL-1)*N
!       *** get indices, density, and attenuation for first layer
        NLYR=NLYRS(ISOL)
        ILYR=1
        M=IXSVP(ILYR,ISOL)

        DEPTH_LOOP: DO I=1,N
!         *** ZI is receiver depth
          ZI=I*DZ
!wyx: begin
          DO WHILE ((ILYR <= NLYR) .AND. (ZI > ZLYR(ILYR,ISOL)))
            ILYR = ILYR + 1
          ENDDO

          IF (ILYR > NLYR) THEN
            ! 接收器位于最底层介质之下: 向外插值

            ILYR = NLYR
            !*** extrapolate
            M = IXSVP(ILYR, ISOL)
            IF (ZSVP(M,ISOL) /= ZSVP(M-1,ISOL)) THEN
              CI = CSVP(M-1,ISOL) + (CSVP(M,ISOL)-CSVP(M-1,ISOL)) *    &
              (ZI-ZSVP(M-1,ISOL)) / (ZSVP(M,ISOL)-ZSVP(M-1,ISOL))
            ELSE
              CI = CSVP(M,ISOL)
            ENDIF
            IF (IEX == 0) THEN
              WRITE(NPU,56)
   56         FORMAT(1X,'WARNING. EXTRAPOLATION OF SVP PERFORMED.')
              IEX = 1
            ENDIF

          ELSE
            ! 接收器位于第 ILYR 层: 内部插值得到ZI深度的当地声速

            IF (ILYR > 1) THEN
              L = IXSVP(ILYR-1, ISOL) + 1
            ELSE
              L = 1
            END IF
            M = IXSVP(ILYR, ISOL)

            !*** determine sound speed CI at depth ZI
            DO J = L, M
              IF (ZI <= ZSVP(J,ISOL)) EXIT
            ENDDO

            !*** interpolate
            CI = CSVP(J-1,ISOL) + (CSVP(J,ISOL)-CSVP(J-1,ISOL)) *    &
            (ZI-ZSVP(J-1,ISOL)) / (ZSVP(J,ISOL)-ZSVP(J-1,ISOL))

          ENDIF

          C1 = CI
          ZG1 = ZLYR(ILYR,ISOL)
          RH1 = RHO(ILYR,ISOL)
          RH1G = RHOG(ILYR,ISOL)
          BETA1 = BETA(ILYR,ISOL)
          BETA1G = BETAG(ILYR,ISOL)

          !*** is receiver on or within 1 dz of interface?
          IF ((ZI + DZ > ZG1) .AND. (ILYR < NLYR)) THEN
            !*** Yes. Get params for medium 2.
            C2 = CSVP(IXSVP(ILYR,ISOL)+1,ISOL)
            ZG2 = ZLYR(ILYR,ISOL)
            RH2 = RHO(ILYR+1,ISOL)
            RH2G = RHOG(ILYR+1,ISOL)
            BETA2 = BETA(ILYR+1,ISOL)
            BETA2G = BETAG(ILYR+1,ISOL)
          ELSE
            !*** Not an interface
            C2 = C1
            ZG2 = ZG1
            RH2 = RH1
            RH2G = RH1G
            BETA2 = BETA1
            BETA2G = BETA1G
          ENDIF
!wyx: end

          RHO1(I,ISOL)=RH1+(ZI-ZG1)*RH1G
          RHO2(I,ISOL)=RH2+(ZI+DZ-ZG2)*RH2G
          BETA1=BETA1+(ZI-ZG1)*BETA1G
          BETA2=BETA2+(ZI+DZ-ZG2)*BETA2G
          IF (BETA1 < 0.0) BETA1=ALPHA*C1/FRQ
          IF (BETA2 < 0.0) BETA2=ALPHA*C2/FRQ
          XN=C0/C1
          XN1(I,ISOL)=CMPLX(XN*XN,XN*XN*BETA1/27.287527)
          XN=C0/C2
          XN2(I,ISOL)=CMPLX(XN*XN,XN*XN*BETA2/27.287527)
        ENDDO DEPTH_LOOP
      ENDDO
END SUBROUTINE
!                                                             **********
!                                                             * PORT2D *
!                                                             **********
      SUBROUTINE PORT2D
!     ******************************************************************
!     * SOLVE FOR FIELD AT PORT WALL BOUNDARY
!     * FIELD AT ADVANCED RANGE IN PORTX
!     ******************************************************************
      USE commons, ONLY: N, AL, AU, AM, DP, PORTX
      IMPLICIT NONE

      CALL TRID( AL(1:N, 1), &
                 AU(1:N, 1), &
                 AM(1:N, 1), &
                 DP(1:N), &
                 PORTX(1:N), N )
!     *** SOLUTION IS IN ARRAY PORTX
      RETURN
      END
!                                                             **********
!                                                             * PORT3D *
!                                                             **********
SUBROUTINE PORT3D
!     ******************************************************************
!     * PORT SIDEWALL BOUNDARY CONDITION
!     ******************************************************************
      USE commons, ONLY: ALFA, DR, DZ, ITYPPW, N, NPU, XK0, RHO1, RHO2, &
                         XN1, XN2, &
                         DP, PORTX, PORTY, SURX, SURY, BOTX, BOTY, U
      IMPLICIT NONE
      REAL    (KIND=8) :: DELTA, R12, RH1, RH2, RSUM, T2
      INTEGER (KIND=4) :: I, IERR
      COMPLEX :: T1,T3,T4,TAL,TAM,TAU,XN12,XNSUM

      SELECT CASE (ITYPPW)
!     *** SET FIELD AT PORT SIDEWALL TO 0.0
      CASE (0)
        DO I=1,N
          PORTX(I)=0.0
          PORTY(I)=0.0
        ENDDO

!     *** ITYPPW = 1 ; USER SUPPLIES SUBROUTINE UPORT3D
      CASE (1)
        CALL UPORT3D

!     *** ITYPPW = 2
!     *** Assumption is that environmental conditions at port sidewall
!     *** boundary are equal to conditions at leftmost solution boundary.
!     *** Solution at boundary at advanced range is solved in 2D.
      CASE (2)
        DELTA=XK0*DR
        T1=CMPLX(.25,.25*DELTA)
        T2=1.0/(XK0*XK0*DZ*DZ)
        T3=T1*T2
        T4=CMPLX(.25,-.25*DELTA)
        RH1=RHO1(1,1)
        RH2=RHO2(1,1)
        R12=1.0/RH1+1.0/RH2
        XN12=T1*((XN1(1,1)-1)/RH1+(XN2(1,1)-1)/RH2)
        RSUM=RH1+RH2
        XNSUM=T1*(XN1(1,1)+XN2(1,1)-2.0)
        TAU=(T3+ALFA+ALFA*T1*(XN2(1,1)-1))/RH2
        TAL=(T3+ALFA+ALFA*T1*(XN1(1,1)-1))/RH1
        TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!       *** COMPUTE RIGHT HAND SIDE D
        DP(1)=TAL*SURY(1)+TAM*PORTY(1)+TAU*PORTY(2)-    &
        (T4*T2+ALFA+ALFA*T4*(XN1(1,1)-1))/RH1*SURX(1)
        DO I=2,N-1
          RH1=RHO1(I,1)
          RH2=RHO2(I,1)
          R12=1.0/RH1+1.0/RH2
          XN12=T1*((XN1(I,1)-1)/RH1+(XN2(I,1)-1)/RH2)
          RSUM=RH1+RH2
          XNSUM=T1*(XN1(I,1)+XN2(I,1)-2.0)
          TAU=(T3+ALFA+ALFA*T1*(XN2(I,1)-1))/RH2
          TAL=(T3+ALFA+ALFA*T1*(XN1(I,1)-1))/RH1
          TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!         *** COMPUTE RIGHT HAND SIDE D
          DP(I)=TAL*PORTY(I-1)+TAM*PORTY(I)+TAU*PORTY(I+1)
        ENDDO
        RH1=RHO1(N,1)
        RH2=RHO2(N,1)
        R12=1.0/RH1+1.0/RH2
        XN12=T1*((XN1(N,1)-1)/RH1+(XN2(N,1)-1)/RH2)
        RSUM=RH1+RH2
        XNSUM=T1*(XN1(N,1)+XN2(N,1)-2.0)
        TAU=(T3+ALFA+ALFA*T1*(XN2(N,1)-1))/RH2
        TAL=(T3+ALFA+ALFA*T1*(XN1(N,1)-1))/RH1
        TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!       *** COMPUTE RIGHT HAND SIDE D
        DP(N)=TAL*PORTY(N-1)+TAM*PORTY(N)+TAU*BOTY(1)-    &
        (T4*T2+ALFA+ALFA*T4*(XN1(N,1)-1))/RH2*BOTX(1)

      CASE DEFAULT
        IERR=1
        WRITE(NPU,50) ITYPPW
   50   FORMAT(1X,'ERROR - PORT SIDEWALL BOUNDARY CONDITION ',I2,    &
        ' NOT IMPLEMENTED')

      END SELECT

      U(1:N, 0) = PORTY(1:N)
END SUBROUTINE
!                                                             **********
!                                                             * PRINTP *
!                                                             **********
      SUBROUTINE PRINTP
!     *****************************************************************
!     *** PRINT PROBLEM PARAMETERS
!     *****************************************************************
!ljh  INCLUDE 'FOR3D.CMN'
      USE commons
      IMPLICIT REAL (KIND=8) (A-H,O-Z)
      IMPLICIT INTEGER (KIND=4) (I-N)

      WRITE(NPU,72)
   72 FORMAT(/1X,'NUSC FOR3D MODEL')
      IF (NDIM == 1) THEN
        WRITE(NPU,*)'2D SOLUTION'
      ELSE
        IF (NDIM == 2) THEN
          WRITE(NPU,*)'N X 2D SOLUTION'
        ELSE
          WRITE(NPU,*)'3D SOLUTION'
        ENDIF
      ENDIF
      IF (ISF == 0) WRITE(NPU,73)
   73 FORMAT(1X,'GAUSSIAN STARTING FIELD')
      IF (ISF == 1) WRITE(NPU,74)
   74 FORMAT(1X,'USER STARTING FIELD')
      IF (ISF == 2) WRITE(NPU,71)
   71 FORMAT(1X,'GREENS STARTING FIELD')
      IF (IHNK /= 0) WRITE(NPU,75)
   75 FORMAT(1X,'STARTING FIELD DIVIDED BY HANKEL FUNCTION')
      IF (ITYPES == 1) WRITE(NPU,76)
   76 FORMAT(1X,'USER SURFACE CONDITION')
      IF (ITYPEB == 1) WRITE(NPU,77)
   77 FORMAT(1X,'USER BOTTOM CONDITION')
      IF (ITYPPW == 1) WRITE(NPU,78)
      IF (NDIM == 3 .AND. ITYPPW == 2)    &
       WRITE(NPU,*)'2D SOLUTION AT PORT SIDEWALL BOUNDARY'
   78 FORMAT(1X,'USER PORT SIDEWALL BOUNDARY CONDITION')
      IF (ITYPSW == 1) WRITE(NPU,79)
      IF (NDIM == 3 .AND. ITYPSW == 2)    &
       WRITE(NPU,*)'2D SOLUTION AT STBD SIDEWALL BOUNDARY'
   79 FORMAT(1X,'USER STARBOARD SIDEWALL BOUNDARY CONDITION')
      IF (NDIM /= 1)ANG=FLDW/NSEC
      WRITE(NPU,80) FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,    &
      ITYPPW,ITYPSW,FLDW,NSEC,ANG,NSOL,NTOT
   80 FORMAT(1X,'FRQ    = ',F8.2,' HZ',/,1X,    &
      'ZS     = ',F8.2,' M',/,1X,    &
      'C0     = ',F8.2,' M/SEC',/,1X,    &
      'ISF    = ',I5,/,1X,    &
      'RA     = ',F8.2,' M',/,1X,    &
      'ZA     = ',F8.2,' M',/,1X,    &
      'N      = ',I5,/,1X,    &
      'IHNK   = ',I5,/,1X,    &
      'ITYPES = ',I5,/,1X,    &
      'ITYPEB = ',I5,/,1X,    &
      'ITYPPW = ',I5,/,1X,    &
      'ITYPSW = ',I5,/,1X,    &
      'FLDW   = ',F8.3,' (WIDTH OF FIELD IN DEGREES)',/,1X,    &
      'NSEC   = ',I5,' SECTORS',/,1X,    &
      'PHI    = ',F8.3,' DEGREES/SECTOR',/,1X,    &
      'NSOL   = ',I5,' SOLUTIONS',/,1X,    &
      'NTOT   = ',I5)
      WRITE(NPU,81) DR,DZ,ANG,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH,RMAX,    &
      DOUGRA,NDIV,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
   81 FORMAT(1X,'DR     = ',F8.2,' M',/,1X,    &
      'DZ     = ',F8.2,' M',/,1X,    &
      'DTH    = ',F8.2,' DEGREES',/,1X,    &
      'WDR    = ',F8.2,' M',/,1X,    &
      'WZ1    = ',F8.2,' M',/,1X,    &
      'WZ2    = ',F8.2,' M',/,1X,    &
      'WDZ    = ',F8.2,' M',/,1X,    &
      'WDTH   = ',F8.2,' DEGREES',/,1X,    &
      'PDR    = ',F8.2,' M',/,1X,    &
      'PDZ    = ',F8.2,' M',/,1X,    &
      'PDTH   = ',F8.2,' DEGREES',/,1X,    &
      'RMAX   =',F9.1,' M',/,1X,    &
      'DOUGRA =',F9.1,' M',/,1X,    &
      'NDIV   = ',I5,/,1X,    &
      'USER 1 = ',F12.4,/,1X,    &
      'USER 2 = ',F12.4,/,1X,    &
      'USER 3 = ',F12.4,/,1X,    &
      'USER 4 = ',F12.4,/,1X,    &
      'USER 5 = ',F12.4,/,1X,    &
      'USER 6 = ',F12.4,/,1X,    &
      'USER 7 = ',F12.4,/,1X,    &
      'USER 8 = ',F12.4,/,1X,    &
      'USER 9 = ',F12.4,/,1X,    &
      'USER 10= ',F12.4,/,1X,    &
      'USER 11= ',F12.4,/,1X,    &
      'USER 12= ',F12.4,/,/)
   84 FORMAT(3X,'LAYER ',I3,3X,'DEPTH = ',F8.1,3X,    &
      'DENSITY = ',F6.2,3X,'GRADIENT = ',F7.4,/,    &
      3X,'BETA = ',F8.3,3X,'GRADIENT = ',F7.4)
!
!     *** PRINT BOTTOM DEPTHS IF REQUESTED
      IF (IBOT == 0) GO TO 86
      DO 85 L=1,NSOL
        TH=THETA(L)*DEG
        IF (NDIM /= 1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC)
!       WRITE(NPU,*)' '
        WRITE(NPU,*)'RELATIVE BEARING =  ',ANG,' DEGREES'
        WRITE(NPU,123) R1(L),Z1(L),R2(L),Z2(L),TH
   85 CONTINUE
   86 IF (ISFLD == 0) GO TO 92
!
!     *** PRINT STARTING FIELD
      WRITE(NPU,87)
   87 FORMAT(/,1X,'STARTING FIELD')
      DO J=1,NSOL
        IF (NDIM /= 1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC)
        WRITE(NPU,*)'RELATIVE BEARING = ',ANG,' DEGREES'
        DO I=1,N
          ZI=I*DZ
          IF (U(I, J) /= 0.0) WRITE(NPU,89) I,ZI,U(I, J)
   89     FORMAT(1X,I4,3X,F10.2,3X,'(',E12.5,2X,E12.5,' )')
        ENDDO
      ENDDO
   92 CONTINUE
!
      IF (ISVP == 0) GO TO 110
!     *** PRINT SVP
      WRITE(NPU,93) RA
   93 FORMAT(//1X,'SOUND VELOCITY PROFILE AT RANGE ',F10.1,' METERS',/)
      DO 102 L=1,NSOL
        NLYR=NLYRS(L)
        IF (NDIM /= 1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC)
!       WRITE(NPU,*)' '
        WRITE(NPU,*)'RELATIVE BEARING = ',ANG,' DEGREES'
!       WRITE(NPU,*)' '
        K=0
        DO 101 ILYR=1,NLYR
          WRITE(NPU,84)ILYR,ZLYR(ILYR,L),RHO(ILYR,L),RHOG(ILYR,L),    &
          BETA(ILYR,L),BETAG(ILYR,L)
!         WRITE(NPU,*)' '
          J=K+1
          K=IXSVP(ILYR,L)
          DO 100 I=J,K
            WRITE(NPU,95) I,ZSVP(I,L),CSVP(I,L)
   95       FORMAT(1X,I4,3X,F8.1,3X,F8.1)
  100     CONTINUE
!         WRITE(NPU,*)' '
  101   CONTINUE
  102 CONTINUE
  110 CONTINUE
      RETURN
  123 FORMAT(1X,'BOTTOM DEPTHS ',/,1X,    &
      'R1    = ',F10.1,' M',/,1X,    &
      'Z1    = ',F8.1,' M',/,1X,    &
      'R2    = ',F10.1,' M',/,1X,    &
      'Z2    = ',F8.1,' M',/,1X,    &
      'THETA = ',F8.1,' DEG')
      END
!                                                             **********
!                                                             *  RHS   *
!                                                             **********
SUBROUTINE RHS
!     *** COMPUTE RIGHT HAND SIDE.
!     *** RA IS NEXT SOLUTION RANGE.
!     *** CAN SAVE TIME BY COMPUTING 2D AND 3D SOLUTIONS IN SEPARATE
!     ***  SUBROUTINES i.e. RHS2 AND RHS3 ETC.
      USE commons, ONLY: N, NDIM, NSOL
      USE commons, ONLY: ALFA, DR, DZ, PHI, RA, XK0, RHO1, RHO2, XN1, XN2, &
                         U, D, SURY, BOTY, PORTY, STBDY
      USE omp_lib
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, L, LM1
      REAL    (KIND=8) :: R12, R34, R56, RH1, RH2, RH3, RH4, RH5, RH6, RSUM, T2
      COMPLEX :: T1,T3,T5,VM,VMM1,VMP1,TBL,TBM,TBU,TBL2,TBM2,TBU2
      COMPLEX :: TAL,TAM,TAU,XN12,XNSUM

!     *** GET SURFACE CONDITION
      CALL SCON3D
!     *** GET BOTTOM CONDITION
      CALL BCON3D
!     *** GET PORT WALL CONDITION
      IF (NDIM == 3) CALL PORT3D
!     *** GET STARBOARD WALL CONDITION
      IF (NDIM == 3) CALL STBD3D

      T1=CMPLX(.25,.25*XK0*DR)
      T2=1.0/(XK0*XK0*DZ*DZ)
      T3=T1*T2
      IF (NDIM == 1 .OR. NDIM == 2 .OR. RA == 0.0) THEN
        T5=0.0
      ELSE
        T5=CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA-.5*DR)*(RA-.5*DR)*PHI*PHI))
      ENDIF
      
     !$OMP PARALLEL DO PRIVATE(LM1,RH3,RH4,R34,TBL,TBU,TBM,RH5,RH6,R56,TBL2)&
	  !$OMP PRIVATE(TBU2,TBM2,VMM1,VM,VMP1,RH1,RH2,R12,XN12,RSUM,XNSUM,TAL,TAU,TAM)
      DO L=1,NSOL
        LM1=MAX(L-1,1)
        DO I=1,N
          RH3=RHO1(I,LM1)
          RH4=RHO1(I,L)
          R34=(1.0/RH3+1.0/RH4)
          TBL=.5*(RH3+RH4)*(T5)/RH3
          TBU=.5*(RH3+RH4)*(T5)/RH4
          TBM=1.0+(-T5)*.5*(RH3+RH4)*R34
          IF (I < N) THEN
            RH5=RHO1(I+1,LM1)
            RH6=RHO1(I+1,L)
            R56=(1.0/RH5+1.0/RH6)
            TBL2=.5*(RH5+RH6)*(T5)/RH5
            TBU2=.5*(RH5+RH6)*(T5)/RH6
            TBM2=1.0+(-T5)*.5*(RH5+RH6)*R56
          ELSE
          ENDIF
          IF (NDIM == 1 .OR. NDIM == 2)  THEN
            TBM=1.0
            TBM2=1.0
            TBL=0.0
            TBU=0.0
            TBL2=0.0
            TBU2=0.0
          ELSE
          ENDIF

          IF (I == 1) THEN
!     *** I=1
            VMM1=TBL*U(I-1,L-1)+TBM*U(I-1,L)+TBU*U(I-1,L+1)
            VM  =TBL*U(I,L-1)+TBM*U(I,L)+TBU*U(I,L+1)
            VMP1=TBL2*U(I+1,L-1)+TBM2*U(I+1,L)+TBU2*U(I+1,L+1)
          ELSEIF (I /= N) THEN
!     *** I=2 TO N-1
            VMM1=VM
            VM  =VMP1
            VMP1=TBL2*U(I+1,L-1)+TBM2*U(I+1,L)+TBU2*U(I+1,L+1)
          ELSE
!     *** I=N
            VMM1=VM
            VM  =VMP1
            VMP1=TBL*U(N+1,L-1)+TBM*U(N+1,L)+TBU*U(N+1,L+1)
          ENDIF

          RH1=RHO1(I,L)
          RH2=RHO2(I,L)
          R12=1.0/RH1+1.0/RH2
          XN12=T1*((XN1(I,L)-1)/RH1+(XN2(I,L)-1)/RH2)
          RSUM=RH1+RH2
          XNSUM=T1*(XN1(I,L)+XN2(I,L)-2.0)
          TAL=(T3+ALFA+ALFA*T1*(XN1(I,L)-1))/RH1
          TAU=(T3+ALFA+ALFA*T1*(XN2(I,L)-1))/RH2
          TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!     *** COMPUTE RIGHT HAND SIDE D
          D(I,L)=TAL*VMM1+TAM*VM+TAU*VMP1
        ENDDO
      ENDDO
END SUBROUTINE
!                                                             **********
!                                                             * SCON3D *
!                                                             **********
SUBROUTINE SCON3D
!     *****************************************************************
!     *** SURFACE BOUNDARY CONDITION SUBROUTINE
!     *****************************************************************
      USE commons, ONLY: ITYPES, NPU, NSOL, &
                         SURX, SURY, U
      IMPLICIT NONE
      INTEGER (KIND=4) :: L, IERR

      SELECT CASE (ITYPES)

!     *** PRESSURE RELEASE SURFACE
      CASE (0)
        DO L=1,NSOL+2
          SURY(L)=0.0
          SURX(L)=0.0
        ENDDO

!     *** USER SURFACE CONDITION
      CASE (1)
        CALL USCON3D

      CASE DEFAULT
        WRITE(NPU,50)ITYPES
   50   FORMAT(1X,'ERROR - SURFACE BOUNDARY CONDITION ',I2,    &
        ' NOT IMPLEMENTED')
        IERR=1

      END SELECT

      U(0, 0:NSOL+1) = SURY(1:NSOL+2)
END SUBROUTINE
!                                                             **********
!                                                             * SFLD3D *
!                                                             **********
SUBROUTINE SFLD3D
!     ***************************************************************
!     *** GENERATE 3D STARTING FIELD
!     ***************************************************************
      USE commons, ONLY: ISF, IHNK, ITYPPW, ITYPSW, N, NSOL, NDIM, NPU
      USE commons, ONLY: C0, U, PORTY, STBDY, DZ, FRQ, PI, RA, XK0, ZS
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, IERR
      REAL    (KIND=8) :: GA, GW, PR, ZM
      COMPLEX          :: HNK,HNKL            ! cljh

      SELECT CASE (ISF)
!     *** ISF=0
!     *** GAUSSIAN STARTING FIELD.
      CASE (0)
        GW=C0/(PI*FRQ)
        GA=1.0/SQRT(GW)
        DO I=1,N
          ZM=I*DZ
          PR=GA*(EXP(-((ZM-ZS)/GW)**2)-EXP(-((-ZM-ZS)/GW)**2))
          U(I, 1:NSOL) = CMPLX(PR,0.0)
        ENDDO
        IF (NDIM /= 1) THEN
          IF (ITYPPW == 2) THEN
            PORTY(1:N)=U(1:N, 1)
          ENDIF
          IF (ITYPSW == 2) THEN
            STBDY(1:N)=U(1:N, 1)
          ENDIF
        ENDIF

!     *** ISF=1
!     *** USER STARTING FIELD
      CASE (1)
      CALL USFLD3D

!     *** ISF=2
!     *** GREENS WIDE ANGLE STARTING FIELD
      CASE (2)
        XK0=2.0*PI*FRQ/C0
        DO I=1,N
          ZM=I*DZ
          PR=SQRT(XK0)*(1.45-0.42*XK0*XK0*(ZM-ZS)**2.0)*    &
          EXP(-(XK0*XK0*(ZM-ZS)**2)/3.0512)
          U(I, 1:NSOL) = CMPLX(PR,0.0)
        ENDDO
        IF (NDIM /= 1) THEN
          IF (ITYPPW == 2) THEN
            PORTY(1:N)=U(1:N, 1)
          ENDIF
          IF (ITYPSW == 2) THEN
            STBDY(1:N)=U(1:N, 1)
          ENDIF
        ENDIF

      CASE DEFAULT
        WRITE(NPU,50)ISF
   50   FORMAT(1X,'ERROR - STARTING FIELD ISF ',I2,' NOT IMPLEMENTED')
        IERR=1

      END SELECT

!     *** DIVIDE STARTING FIELD BY HANKEL FUNCTION IF REQUESTED BY USER
      IF (IHNK /= 0) THEN
        HNK=HNKL(XK0*RA)
        U(1:N, 1:NSOL) = U(1:N, 1:NSOL) / HNK
        PORTY(1:N)=PORTY(1:N)/HNK
        STBDY(1:N)=STBDY(1:N)/HNK
      ENDIF
END SUBROUTINE
!                                                             **********
!                                                             * STBD2D *
!                                                             **********
SUBROUTINE STBD2D
!     ******************************************************************
!     * SOLVE FOR FIELD AT STBD WALL BOUNDARY
!     * FIELD AT ADVANCED RANGE IN STBDX
!     ******************************************************************
      USE commons, ONLY: N, NSOL, AL, AU, AM, DS, STBDX
      IMPLICIT NONE

      CALL TRID( AL(1:N, NSOL), &
                 AU(1:N, NSOL), &
                 AM(1:N, NSOL), &
                 DS(1:N), &
                 STBDX(1:N), N )
!     *** SOLUTION IS IN ARRAY STBDX
END SUBROUTINE
!                                                             **********
!                                                             * STBD3D *
!                                                             **********
SUBROUTINE STBD3D
!     ******************************************************************
!     * STARBOARD SIDEWALL BOUNDARY CONDITION
!     ******************************************************************
      USE commons, ONLY: N, ITYPSW, NPU, NSOL
      USE commons, ONLY: DR, XK0, DZ, DS, ALFA, RHO1, RHO2, XN1, XN2, &
                         SURX, SURY, STBDX, STBDY, BOTX, BOTY, U
      IMPLICIT NONE

      REAL    (KIND=8) :: DELTA, R12, RH1, RH2, RSUM, T2
      INTEGER (KIND=4) :: I, IERR, M
      COMPLEX          :: T1,T3,T4,TAL,TAU,TAM,XN12,XNSUM

      SELECT CASE (ITYPSW)
!     *** ITYPSW = 0
!     *** SET FIELD AT STBD SIDEWALL TO 0.0
      CASE (0)
        DO I=1,N
          STBDX(I)=0.0
          STBDY(I)=0.0
        ENDDO

!     *** ITYPSW = 1 ; USER SUPPLIES SUBROUTINE USTBD3D
      CASE (1)
        CALL USTBD3D

!     *** ITYPSW = 2
!     *** ASSUMPTION IS THAT ENVIRONMENTAL CONDITIONS AT STBD SIDEWALL
!     *** BOUNDARY ARE EQUAL TO CONDITIONS AT RIGHTMOST SOLUTION BOUNDARY.
!     *** SOLUTION AT BOUNDARY AT ADVANCED RANGE IS SOLVED IN 2D.
      CASE (2)
        DELTA=XK0*DR
        T1=CMPLX(.25,.25*DELTA)
        T2=1.0/(XK0*XK0*DZ*DZ)
        T3=T1*T2
        T4=CMPLX(.25,-.25*DELTA)
        M=(NSOL-1)*N
        RH1=RHO1(1,NSOL)
        RH2=RHO2(1,NSOL)
        R12=1.0/RH1+1.0/RH2
        XN12=T1*((XN1(1,NSOL)-1)/RH1+(XN2(1,NSOL)-1)/RH2)
        RSUM=RH1+RH2
        XNSUM=T1*(XN1(1,NSOL)+XN2(1,NSOL)-2.0)
        TAU=(T3+ALFA+ALFA*T1*(XN2(1,NSOL)-1))/RH2
        TAL=(T3+ALFA+ALFA*T1*(XN1(1,NSOL)-1))/RH1
        TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!       *** COMPUTE RIGHT HAND SIDE D
        DS(1)=TAL*SURY(NSOL+2)+TAM*STBDY(1)+TAU*STBDY(2)-    &
        (T4*T2+ALFA+ALFA*T4*(XN1(1,NSOL)-1))/RH1*SURX(NSOL+2)
        DO I=2,N-1
          RH1=RHO1(I,NSOL)
          RH2=RHO2(I,NSOL)
          R12=1.0/RH1+1.0/RH2
          XN12=T1*((XN1(I,NSOL)-1)/RH1+(XN2(I,NSOL)-1)/RH2)
          RSUM=RH1+RH2
          XNSUM=T1*(XN1(I,NSOL)+XN2(I,NSOL)-2.0)
          TAU=(T3+ALFA+ALFA*T1*(XN2(I,NSOL)-1))/RH2
          TAL=(T3+ALFA+ALFA*T1*(XN1(I,NSOL)-1))/RH1
          TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!         *** COMPUTE RIGHT HAND SIDE D
          DS(I)=TAL*STBDY(I-1)+TAM*STBDY(I)+TAU*STBDY(I+1)
        ENDDO
        RH1=RHO1(N,NSOL)
        RH2=RHO2(N,NSOL)
        R12=1.0/RH1+1.0/RH2
        XN12=T1*((XN1(N,NSOL)-1)/RH1+(XN2(N,NSOL)-1)/RH2)
        RSUM=RH1+RH2
        XNSUM=T1*(XN1(N,NSOL)+XN2(N,NSOL)-2.0)
        TAU=(T3+ALFA+ALFA*T1*(XN2(N,NSOL)-1))/RH2
        TAL=(T3+ALFA+ALFA*T1*(XN1(N,NSOL)-1))/RH1
        TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12
!       *** COMPUTE RIGHT HAND SIDE D
        DS(N)=TAL*STBDY(N-1)+TAM*STBDY(N)+TAU*BOTY(NSOL+2)-    &
        (T4*T2+ALFA+ALFA*T4*(XN1(N,NSOL)-1))/RH2*BOTX(NSOL+2)

      CASE DEFAULT
        IERR=1
        WRITE(NPU,50) ITYPSW
   50   FORMAT(1X,'ERROR - STARBOARD SIDEWALL BOUNDARY CONDITION ',I2,    &
        ' NOT IMPLEMENTED')

      END SELECT

      U(1:N, NSOL+1) = STBDY(1:N)
END SUBROUTINE
!                                                             **********
!                                                             * SVP3D  *
!                                                             **********
SUBROUTINE SVP3D
!     ******************************************************************
!     *** SOUND VELOCITY PROFILE SUBROUTINE
!     ******************************************************************
!     *** CALLING PROGRAM SUPPLIES: NOTHING
!     *** SVP SUBROUTINE RETURNS:
!         NLYRS - NUMBER OF LAYERS. LAYER 1 IS WATER. OTHERS ARE SEDIMENT
!         ZLYR - ARRAY - DEPTH OF EACH LAYER. FIRST IS DEPTH OF WATER.
!         RHO  - ARRAY - DENSITY OF EACH LAYER. GRAMS/CUBIC CM
!         RHOG - ARRAY - DENSITY GRADIENT.
!         BETA - ARRAY - ATTENUATION IN EACH LAYER. DB/WAVELENGTH
!         BETAG- ARRAY - ATTENUATION GRADIENT.
!         IXSVP- ARRAY - CONTAINS POINTERS. POINTS TO LAST VALUE OF SVP
!                IN CORRESPONDING LAYER. SVP IS STORED IN ARRAYS ZSVP
!                AND CSVP. IXSVP(1,L) POINTS TO LAST SVP POINT IN WATER
!                ALONG SECTOR BOUNDARY L. L=1,NSOL.
!                IXSVP(NLYR,L) POINTS TO LAST SVP POINT IN BOTTOM-MOST LAYER.
!         NSVP - NUMBER OF POINTS IN ZSVP AND CSVP. ZSVP AND CSVP
!                CONTAIN THE PROFILES FOR ALL LAYERS.
!         ZSVP - ARRAY - SVP DEPTHS - METERS
!         CSVP - ARRAY - SOUND SPEED - METERS/SEC
!     ******************************************************************
!
      USE commons, ONLY: NLYRS, IXSVP, NIU, NLYR, NSOL, NSVP
      USE commons, ONLY: ZLYR, RHO, RHOG, BETA, BETAG, ZSVP, CSVP
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, J, L, LASTL, LNLYR
      REAL    (KIND=8) :: CV, ZV

!     *** READ NUMBER OF LAYERS
      DO L=1,NSOL
        NSVP=0
        READ(NIU,*,END=100) NLYRS(L)
        NLYR=NLYRS(L)
!
!     *** IF NLYR IS NEGATIVE, REPEAT LAST PROFILE UNTIL NSOL PROFILES
!     *** HAVE BEEN ENTERED.
        IF (NLYR < 0) GO TO 200
        LNLYR=NLYR
        LASTL=L
!     *** FIRST LAYER IS WATER. OTHERS ARE SEDIMENT.
        DO 55 I=1,NLYR
!     *** READ DEPTH OF LAYER, DENSITY AND ATTENUATION.
          READ(NIU,*,END=100) ZLYR(I,L),RHO(I,L),RHOG(I,L),BETA(I,L),    &
          BETAG(I,L)
!         *** READ PROFILE.
   50     READ(NIU,*,END=100) ZV,CV
          NSVP=NSVP+1
          ZSVP(NSVP,L)=ZV
          CSVP(NSVP,L)=CV
          IF (ZV < ZLYR(I,L)) GO TO 50
          IF (ZV > ZLYR(I,L)) GO TO 100
          IXSVP(I,L)=NSVP
   55   CONTINUE
      ENDDO
      RETURN
!
!     *** ERROR EXIT
  100 NSVP=0
      RETURN
!
!     *** REPEAT LAST PROFILE UNTIL NSOL PROFILES ENTERED.
  200 CONTINUE
      NLYR=LNLYR
      NSVP=IXSVP(NLYR,LASTL)
      DO L=LASTL+1,NSOL
        NLYRS(L)=NLYR
        DO I=1,NLYR
          ZLYR(I,L)=ZLYR(I,LASTL)
          RHO(I,L)=RHO(I,LASTL)
          RHOG(I,L)=RHOG(I,LASTL)
          BETA(I,L)=BETA(I,LASTL)
          BETAG(I,L)=BETAG(I,LASTL)
          IXSVP(I,L)=IXSVP(I,LASTL)
        ENDDO
        DO J=1,NSVP
          ZSVP(J,L)=ZSVP(J,LASTL)
          CSVP(J,L)=CSVP(J,LASTL)
        ENDDO
      ENDDO
END SUBROUTINE

!                                                             **********
!                                                             * TRID   *
!                                                             **********
SUBROUTINE TRID(XL, XU, XM, B, X, N)
      IMPLICIT NONE
      INTEGER, INTENT( IN ) :: N
      COMPLEX, INTENT( IN ) :: XL(1:N), XU(1:N), XM(1:N), B(1:N)
      COMPLEX, INTENT( OUT) :: X(1:N)
      COMPLEX :: BTA(1:N), GAMMA(1:N)
      INTEGER :: I

      BTA(1) = XM(1)
      DO I = 2, N
        BTA(I) = XM(I) - XL(I) * XU(I-1) / BTA(I-1)
      ENDDO
      DO I = 1, N
        BTA(I) = 1.0 / BTA(I)
      ENDDO

      GAMMA(1) = B(1) * BTA(1)
      DO I = 2, N
        GAMMA(I) = ( B(I) - XL(I) * GAMMA(I-1) ) * BTA(I)
      ENDDO

      X(N) = GAMMA(N)
      DO I = N, 2, -1
        X(I-1) = GAMMA(I-1) - XU(I-1) * X(I) * BTA(I-1)
      ENDDO
END SUBROUTINE


!                                                             **********
!                                                             * TWOSTEP*
!                                                             **********
SUBROUTINE TWOSTEP
!     ******************************************************************
!     * SOLVE V
!     * SOLVE U
!     * FIELD AT ADVANCED RANGE IN U
!     ******************************************************************
      USE commons, ONLY: ITYPPW, ITYPSW, N, NDIM, NSOL
      USE commons, ONLY: AL, AU, AM, BL, BU, BM, D, U, &
                         PORTX, PORTY, STBDX, STBDY, &
                         ALFA, DR, DZ, PHI, RA, XK0, &
                         RHO1, RHO2, SURX, BOTX, XN1, XN2
      USE omp_lib
      IMPLICIT NONE
      INTEGER (KIND=4) :: I, L, LWYX
      REAL    (KIND=8) :: R34, RH3, RH4
      COMPLEX          :: T1,T2,T3,T4,T5,TBL,TBM,TBU,TAL,TAU
      COMPLEX          :: VBOTX,VSURX

      IF (NDIM == 3 .AND. ITYPPW == 2) CALL PORT2D
      IF (NDIM == 3 .AND. ITYPSW == 2) CALL STBD2D

      T1=CMPLX(.25,.25*XK0*DR)
      T2=CMPLX(.25,-.25*XK0*DR)
      T3=T2/(XK0*XK0*DZ*DZ)
      T4=T1/(XK0*XK0*DZ*DZ)
      IF (NDIM == 1 .OR. NDIM == 2 .OR. RA == 0.0) THEN
        T5=0.0
      ELSE
        T5=-CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA-.5*DR)*(RA-.5*DR)*PHI*PHI))
      ENDIF
      
       !$OMP PARALLEL DO PRIVATE(LWYX,TBL,TBU,TBM,RH3,RH4,R34,VSURX,TAL,VBOTX,TAU)
      DO L=1,NSOL
        LWYX=MAX(L-1,1)

        IF (NDIM == 1 .OR. NDIM == 2) THEN
          TBL=0.0
          TBU=0.0
          TBM=1.0
        ELSE
          RH3=RHO1(1,LWYX)
          RH4=RHO1(1,L)
          R34=(1.0/RH3+1.0/RH4)
          TBL=.5*(RH3+RH4)*(T5)/RH3
          TBU=.5*(RH3+RH4)*(T5)/RH4
          TBM=1.0+(-T5)*.5*(RH3+RH4)*R34
        ENDIF
        VSURX=TBL*SURX(L)+TBM*SURX(L+1)+TBU*SURX(L+2)
        TAL=(T3+ALFA+ALFA*T2*(XN1(1,L)-1))/RHO1(1,L)
        D(1,L)=D(1,L)-TAL*VSURX

        IF (NDIM == 1 .OR. NDIM == 2) THEN
          TBL=0.0
          TBU=0.0
          TBM=1.0
        ELSE
          RH3=RHO1(N,LWYX)
          RH4=RHO1(N,L)
          R34=(1.0/RH3+1.0/RH4)
          TBL=.5*(RH3+RH4)*(T5)/RH3
          TBU=.5*(RH3+RH4)*(T5)/RH4
          TBM=1.0+(-T5)*.5*(RH3+RH4)*R34
        ENDIF
        VBOTX=TBL*BOTX(L)+TBM*BOTX(L+1)+TBU*BOTX(L+2)
        TAU=(T3+ALFA+ALFA*T2*(XN2(N,L)-1))/RHO2(N,L)
        D(N,L)=D(N,L)-TAU*VBOTX

        CALL TRID( AL(1:N, L), &
                   AU(1:N, L), &
                   AM(1:N, L), &
                   D (1:N, L), &
                   D (1:N, L), N )
      ENDDO
!     *** SOLUTION IS IN ARRAY U
      IF (NDIM == 1 .OR. NDIM == 2) RETURN

    !$OMP PARALLEL DO
      DO I = 1, N
        D(I, 1)    = D(I, 1) - PORTX(I) * T5
        D(I, NSOL) = D(I, NSOL) - STBDX(I) * T5

        CALL TRID( BL(I, 1:NSOL), &
                   BU(I, 1:NSOL), &
                   BM(I, 1:NSOL), &
                   D (I, 1:NSOL), &
                   U (I, 1:NSOL), NSOL )
      ENDDO
!     *** SOLUTION IN V BUT MUST BE REORDERED.
!     *** REORDER SOLUTION
      IF (ITYPPW == 2) PORTY(1:N)=PORTX(1:N)
      IF (ITYPSW == 2) STBDY(1:N)=STBDX(1:N)
!     *** SOLUTION FIELD AT ADVANCED RANGE IS IN ARRAY U.
      RETURN
END SUBROUTINE
